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ABSTRACT 



The problem of roll, sway and yaw motions of surface ships is considered. A 
mathematical model is developed which consists of the nonlinear maneuvering equations 
and incorporates cross coupling between sway force, yaw moment and the roll angle induced 
during a steady turn. The hydrodynamic derivatives and coefficients of a typical container 
ship were used as the base-line study model. The coupled system of nonlinear algebraic 
equations is formulated and solved to predict the steady state roll angle, sway velocity and 
turning rate as a function of the rudder angle. The results are then compared to that of the 
decoupled systems currently employed. A local perturbation is implemented in the vicinity 
of the above steady states to investigate dynamic stability of motion. Sensitivity analysis with 
respect to important design parameters such as speed loss during turning, approach speed, 
transverse metacentric height and trim is performed. Results demonstrate the significance 
of the coupling between roll, sway and yaw and the need to incorporate similar studies in 
the ship design and analysis process. 
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I. INTRODUCTION 



A. GENERAL 

Any study of the motion of a ship through water requires the use of models. 
The model could be a mathematical representation which attempts to describe the 
situation symbolically, searching for a quantitative solution, or a physical model 
which strives to replicate or simulate the actual conditions for the purpose of 
collecting empirical data. The two can be combined such that they complement each 
other, or most often, physical models are used to test the validity and accuracy of 
mathematical models. However, both models fail to completely describe the physical 
phenomenon. The physical model falls short by the inability to simultaneously satisfy 
all the parameters describing the problem. The mathematical model suffers by 
departures from the true state through assumptions made during the development 
of the equations [Ref. 1]. Notwithstanding these shortcomings, the use of 
models in predicting a ships response to its environment, i.e., wind and waves, and 
its own propulsive and control systems is essential in the design process to ensure 
safety of the ship, performance of its mission, its survivability in extreme conditions 
and its efficiency during normal operating conditions. [Ref. 1] 

In general, the examination of ship motions is separated into: 1) steering and 
maneuverability (calm water), and 2) seakeeping (waves, current), with motion 
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stability (absence of external excitation) and control (external excitation) some of the 
primary concerns. [Ref. 1] 

The basis of the mathematical model is rigid body dynamics, hence Newton’s 
laws of motion: [Ref. 2] 

F - — x (linear momentum ) 

* j (11) 

M = — x (angular momentum) 

The ship is free to rotate (roll, pitch, yaw) and translate (surge, sway, heave) in all 
six degrees of freedom and the forces and moments acting on the vehicle are 
comprised of: [Ref. 2] 

1. Hydrodynamic forces and moments on the bare hull, appendages, rudder and 
propeller. 

2. Inertial reaction forces and moments. 

3. Wind, waves and currents. 

4. External forces. 

Moreover, these forces and moments are dependent upon properties of the rigid 
body (e.g., geometry, mass, center of gravity), its orientation (i.e., to inertial frame 
of reference and to the body fixed frame of reference), its dynamics (velocities, 
accelerations, propeller speed, rudder deflection and deflection rate), and the fluid 
properties (e.g., density, viscosity, pressure, energy). [Ref. 1] 
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At this point, the mathematical model has rapidly developed into a very 
complex and cumbersome system and requires simplification in order to be of 
practical use. Choice of the body fixed coordinate system to take advantage of any 
symmetries is one simplification used almost exclusively. Eda, et al, [Ref. 3] 
found in one study of bulbous bow type ships that the slight asymmetry of the 
underwater hull form increased significantly in roll such that this simplification was 
not valid. Another reduction in complexity often used is to separate the lateral 
motions (sway, yaw, roll) from the longitudinal motions (surge, pitch, heave) 
[Refs. 4, 5, 6]. Later studies [Refs. 3, 7, 8, 9, 10] account for the loss 
of speed during a turn by including the surge equation with the lateral motions for 
maneuvering studies. Seakeeping studies [Refs. 8, 11, 12, 13] heretofore have 
separated out the roll equations using a one degree of freedom system only, coupling 
effects were neglected as very small or incorporated as forcing terms. Rutgersson and 
Ottosson [Ref. 8] superimposed the maneuvering model motions with the seakeeping 
motions. 

As can be ascertained from this discussion, the forces, moments, velocities and 
accelerations acting upon the vessel are highly nonlinear, involving complex coupling 
of terms. Linear theory has been successful in predicting and analyzing directionally 
stable ships with controls fixed for small perturbations [Ref. 14]. In addition, the 
linear hydrodynamic forces and moments have been reduced to semi-empirical 
equations based on ship design characteristics such as block coefficient, length, 
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breadth and draft [Refs. 2, 15]. However, there has been no completely analytical 
procedure which predicts the nonlinear hydrodynamic forces and moments [Ref. 2]. 

It is necessary therefore, to combine the mathematical model with the physical 
model such that the physical model is used to measure these nonlinear forces and 
moments and incorporate the results in the mathematical model. "The functions that 
may approximate the hydrodynamic forces and moments can be expressed formally 
by Taylor expansion around the state of equilibrium with respect to the quantities 
affecting the different motions, such as the axial speed u, turning rate r, rudder angle 
8, etc. The polynomial coefficients thus describing the hydrodynamic forces and 
motions can be determined from captive model tests at different kinds of basic 
motions [Ref. 8]. " A complete treatment of the captive model test procedure can be 
found in Principals of Naval Architecture, Vol. Ill [Ref. 2], and will not be discussed 
here. The Taylor expansion and subsequent model development is elaborated in the 
following Chapters. 

Captive model tests for the lateral motions are more difficult than for the 
longitudinal motions due to the large sizes required to adequately represent the 
physical model. Consequently, there exists little data for the third order coefficients 
necessary to produce semi-empirical equations for the nonlinear hydrodynamic forces 
and moments arising from these lateral motions similar to those which exist for the 
first order terms. One study by Inoue, et al, [Ref. 15] has produced formulations for 
limited nonlinear terms in yaw and sway based on drift angle and turning rate. The 
yaw moment coefficients so produced are not in good agreement with the model test 
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data under all conditions of loading and motions induced. Another study done by 
Son and Nomoto [Ref. 9], utilizes model test data for the third order hydrodynamic 
coefficients of a SR- 108 container ship. In the absence of such stmi-empirical 
formulations for the forces and moments in the equations of motion and access to 
model test facilities, the aforementioned data, together with vehicle design 
characteristics, can be used in a mathematical model to ascertain particular 
information about the coupling of yaw, sway and roll motions. 

B. SCOPE OF THIS THESIS 

This research is intended as a bridge between the static roll restoring moment, 
controls fixed analysis of the coupled lateral motion problem and a fully coupled 
treatment of a ships roll response to maneuvering in a seaway. As discussed 
previously, current studies in the area of maneuvering and steering range from three 
degree of freedom models, to four degree of freedom models which account for the 
speed loss during a turn. These studies, however, use the equilibrium point in the 
Taylor expansion of the representative hydrodynamic forces and moments, as straight 
upright, zero rudder deflection, such that v 0 =r 0 =</> 0 =S 0 =0 and all perturbations and 
nonlinearity effects are evaluated based on this single nominal point. This research 
has as its basic premise that the nominal point for the Taylor expansion changes 
significantly with rudder angle. As a result of this change, the stability characteristics 
of the ship may also change appreciably, a fact which needs to be established for 
successful prediction of ship maneuvering response in a seaway. 
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Therefore, the objective is to determine the significance of the coupling of roll 
into sway and yaw motions as a function of rudder angle. In analyzing this coupled 
problem, the effects of roll into sway/yaw and conversely, sway/yaw into roll were 
investigated. 

Development of the mathematical model for the coupled yaw, sway, and roll 
equations of motion along with models for the rudder forces and moments, the 
righting moment, and the hydrodynamic forces and moments is presented in Chapter 
II. The final form of the model being reduced to a set of four nonlinear, coupled, 
first order differential equations. The additional equation is derived from kinematics 
to relate the roll angle and the roll rate. 

In order to utilize Taylor expansion for formulation of the functions which 
approximate the hydrodynamic forces and moments, the state of equilibrium for the 
motion parameters, v, r, and </>, is established. This is accomplished through steady 
state analysis. In Chapter III, this analysis is first done for the linear case to obtain 
an initial approximation for the case of near zero rudder deflection (S=0). This 
linear steady state solution is incorporated into the nonlinear steady state analysis 
in an iterative process which produced the equilibrium values for yaw rate, sway 
velocity, and roll angle as a function of the rudder deflection in a steady turn. As the 
rudder angle is incrementally increased, the steady state values from the previous 
rudder angle are used as the initial approximations to initiate the iterations for the 
current rudder deflection. 
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Once the steady state solutions have been computed, their stability properties 
are established by local perturbation of the coupled equations of motion in the 
neighborhood of each nominal point. This procedure, described in Chapter IV, yields 
a generalized eigenvalue problem. The solutions to the eigenvalue problem provide 
information as to the stability properties of the nominal points and thus characterizes 
the stability of the coupled sway, yaw and roll motions. Demonstration of the 
coupling effect on the lateral motions, is shown utilizing comparison to two 
independent systems: 

1. A sway/yaw model represented by two nonlinear equations decoupled from 
roll in that the hydrodynamic forces and moments are functions of sway 
velocity (v) and yaw rate (r) only. 

2. A roll model represented by the nonlinear roll equation and the kinematic 
roll rate equation decoupled from sway and yaw in that the hydrodynamic 
forces and moments are functions of roll angle (<f>) and roll rate (p) only. 

The directional and roll stability characteristics in a steady turn are established 
by the degree to which the real part of the complex eigenvalues is negative and the 
magnitude of the roll damping ratio. 

Chapter V relates the sensitivity and stability of design parameters such as 
metacentric height and trim to the coupling between the lateral motions. Results are 
presented in terms of rudder angle with a qualitative summary at the end of this 
chapter. 

Finally, Chapter VI summarizes the conclusions and presents recommendations 
for further research and utilization of the mathematical model herein developed. 
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II. MATHEMATICAL MODEL 



A. EQUATIONS OF MOTION 

For ship motions, the full nonlinear equations of motion in a body fixed axis 
system are shown in Equations (2.1). These equations reflect an assumption of 
symmetry along the longitudinal axis (i.e., y G =0). The orthogonal, right-hand axes for 
the body fixed system is shown in Figure 2.1. The parameters are defined in Table 
I. 



Surge: 


m [ii -x G (q 2 +r 2 ) +z G (q +pr) -rv+wq] = X 


(2.1a) 


Sway: 


m[v+x G (pq+f)+z G (.qr-p)+ur-wp] = Y 


(2.1b) 


Heave: 


m[w+x G (rp-q)-z G (p 2 +q 2 )+pv-qu] = Z 


(2.1c) 


„ I ,P + ( I r I ^ r ~ I xy(M + ^ ~mz G (v+ur~wp) 
Roll: 

= K-AGZ(<k) 


(2- Id) 


Pitch: 


I y q -(I z ~I x )rp -IJr 2 -p 2 ) +m[z G (u +qw-rv) - 

x^w+pv-qu)] =M 


(2T e) 


Yaw: 


Jj+Vy-Qpq + I xz (qr-p) +mx G (v+ur-wp) = N 


(2- If) 
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B. SIMPLIFYING ASSUMPTIONS 

For the purpose of this investigation, a mathematical model having three 
degrees of freedom vice the six degrees of freedom depicted in Equations (2.1) was 
used. The following simplifying assumptions were made: 



1. The rotational velocity and acceleration about the y-axis are zero. ( q = 0 
and q = 0 ) 

2. The translational velocity and acceleration in the z direction are zero. ( w 
= 0 and w = 0 ) 

3. The vertical heave and pitch motions are decoupled from the horizontal plane 
motions. 

4. The product of inertia 1^ is very small and can be neglected. 

5. The surge equation is substituted by an algebraic equation which is a function 
of u, V, and 5. 

6. The longitudinal center of gravity, (LCG) and the longitudinal center of 
buoyancy, (LCB) are at midship. 
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7. The vertical center of gravity, (VCG) is on the centerline. 



8. The only important forces and moments acting on the ship induced by the 
rudder are those due to rudder deflection. Forces and moments due to 8 and 
8 are negligible. [Ref. 2] 

TABLE I. EQUATIONS OF MOTION PARAMETERS 



Parameter 


Description 


x,y,z 


Distance along the principal axes 


u,v,w 


Translational velocity components of ship relative to fluid 
along body axes 


p.q.r 


Rotational velocity components of ship relative to inertial 
reference system along body axes 


X,Y,Z 


Hydrodynamic force components along body axes 


K,M,N 


Hydrodynamic moment components along body axes 




xp yaw angle: bow to starboard positive 

6 pitch angle: bow up positive 

<f> roll angle: starboard down positive 


m 


Mass of ship 


Xo.yc.Zo 


Coordinates of the center of gravity in the body axis 
system 


Ix.Iy.Iz 


Moments of inertia about the body axis system 


Ixz.Iyz.Ixy 


Products of inertia about the body axis system 


L 


Displacement volume of ship 


A 


Displacement weight of ship 


GZ(<f>) 


Righting moment as a function of roll angle 


n 


Rudder angle in radians 


V 


Initial velocity of ship 


p 


Mass density of sea water 


L 


Ship length between perpendiculars (LBP) 
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Applying these assumptions to Equations (2.1), the three degree of freedom 



equations of motion for the model used are given in Equations (2.2). 



Sway: 


miv+xj-ZfP+u^ = Y 


(2.2a) 


Roll: IJ>- 


-mz G (v+ur) = AT- AGZ(4>) 


(2.2b) 


Yaw: 


/jf+mx^v+ur) = N 


(2.2c) 



Equations (2.2) are nondimensionalized for ease of working between model test 
data and actual ship test data using the relationships shown in Table II. 

To demonstrate the use of these nondimensionalizing terms, substitution of the 
appropriate values from Table II into Equation (2.2a) yields: 
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TABLE II. NONDIMENSIONAL PARAMETER RELATIONSHIPS 



/ V 

V = — 

V 


• / • 
V = V 


f L \ 




y' = * 


' 1 ) 




{ v 2 i 


Ipi’K 2 

12 J 








** 

ii 




'o' 




N' = N 


1 ' 






Ipl 3 *’ 2 

U , 


4 > 7 = 4 > 


4,' = < 


(t) 




K' = K 


1 ' 
Ipt»F 2 

U , 




^3 

II 

^ 1 


p 1 = p 






/ “ 
w = — 

V 


m' = m 


1 ' 

| Pi ’ 
V z 




X G = X 


« 




M 

O 

II 

1 ^ 








X*" 

II 

w*- 


r 1 \ 
i pi5 




h - V 

< 


1 ] 
I PL ’ 

z y 


AGZ( 4 >)' 

/ 

AGZ( 4 >) 

* 


1 ] 
IpI’K 2 

,2 P J 
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simplifying by factoring out (V 2 /L): 



— dL 2 V*\ m' \v , +x!~j , -zij>'+u , r f ] = (•^pI 2 K 2 j Y' 




and Equation (2.2a) becomes Equation (2.3) in nondimensional form. 



m' [ + x' G r' - z' G p' + u'r' ] = Y 1 



(2.3) 



Since Equation (2.3) and Equation (2.2a) are of the same form, the prime notation 
will be dropped and all equations will be considered as represented in 
nondimensional form unless otherwise indicated. 

C. FORCE AND MOMENT REPRESENTATION 

Using the method of Abkowitz and Strom-Tejsen [Ref. 2], the sway force, roll 
moment and yaw moment can be expressed as Equations (2.4). 



A third order Taylor expansion of f u f 2 , and f 3 can be expressed as Equations 



Sway Force : Y = / 1 (u,v,r,ii,v,r,<J>,<}>,fl) (2.4a) 

Roll Moment. K = f 2 (u,v,r,u,v,r,<}),<j>,4>,&) (2.4b) 

Yaw Moment : N = f 2 (u,v,r,ti,v,r,4>,4>,6) (2.4c) 



(2.5). 
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(2.5b) 



(2.5c) 



f 2 = K v v +AT vvv v 3 2 r +K/ 3 'K^vr 2 +K^ +K$ +K $ + 

tf^v 2 4> + ^W> 2 +**/4> 2 +*(6) +* 0 - A GZ (4») 

/ 3 = N v v+N^y+N vvv v 3 +N vn yr 1 +N r r+N/+N rn r 2 +N vw v 2 r+ 

+/tyj> + ^v 2 <{, ^v^+N^r 2 * +Af nM / , 4> 2 +N(6) +N 0 

The coefficients in Equations (2.5) are the hydrodynamic coefficients and 

represent the partial derivative with respect to the subscripted variable. For example, 

'&y) 



v dY , v a 

Y v means — and Y vrr means — 
dv dv 






The terms Y 0 , Kq, and N 0 are the sway force, 
roll moment and yaw moment induced by the propeller respectively. In a similar 
manner Y(S), K(S) and N(S) represent the force and moments induced by the 
rudder. The righting moment due to the static stability of the ship is represented by 
the AGZ ((f>) term. 



D. RUDDER FORCE AND MOMENT REPRESENTATION 

The expressions used to determine the rudder force and moments were taken 
from Son and Nomoto [Ref. 9] and are presented in Equations (2.6) and (2.7). The 
parameters are defined in Table III. 



m = 


- (1 + a u )F N cos(6) 


(2.6a) 


m) = 


(i + a u )z R F N c os(6) 


(2.6b) 


W) - 


- ( x R + a^Fff cos(6) 


(2.6c) 
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TABLE III. RUDDER FORCE AND MOMENT PARAMETERS 



Parameter 


Description 


a H 


Rudder to hull interaction coefficient 


F N 


Normal force action on the rudder 


Z R 


z coordinate of point on which rudder force Y 6 acts 


Xr 


x coordinate of point on which rudder force Y s acts 


X H 


x coordinate of point on which normal force F N acts 


A 


Rudder aspect ratio 


V R 


Rudder area 


Vr 


Effective rudder inflow velocity 


a M 


Effective rudder inflow angle 


UrAr 


Components of rudder effective inflow velocity 


e 


constant in Equation (2.7c) 


u p 


Effective propeller inflow velocity 


m 


constant in Equation (2.7c) 


K T 


Thrust coefficient 


J 


Advance coefficient 


n 


Number of propeller revolutions per second 


D 


Propeller diameter 


w p 


Effective propeller wake fraction 


J 


constant in Equation (2.7e) 


*P 


x coordinate of propeller position 


^pv>^pr 


propeller flow rectification coefficients 


Y 


flow rectification coefficients 


^5r’^5rrr’^5rrv 


rudder wake coefficients 
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V r sinCoc^) 



(2.7a) 




(2.7b) 



u 



R 




8 kK T 
it J 2 



(2.7c) 



uV 
J = - £ — 




(2.7d) 


nD 

u p = u [(1 -wp +x{(v+x p r) 2 +c /TV ,v+c pr r}] 


(2.7e) 


V A = Y V+C 4/ +C 8m /,3+C 8rr/ 2 V 


(2.7f) 


ctjj = d + tan -1 




(2-7g) 



E. RIGHTING MOMENT REPRESENTATION 

For the initial part of the righting arm curve and for wall sided ships, Equation 
(2.8) gives a good approximation for the function GZ (<£). [Ref. 16] 

GZ(<j>) = GMsin(4>) + — tan 2 (<j>)sin(<}>) (2.8) 

2 

GM = the transverse metacentric height 
where BM = the transverse metacentric radius 

<J )= the roll angle 



Equation (2.8) expresses the fact that for most ships of fairly rectangular 
midship section, the GZ(<£) curve exhibits the typical characteristics of a hardening 
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spring restoring moment. The spring constant is an increasing function of the roll 
angle (f>. To see this, we use the Taylor series expansion up to third order for the 
trigonometric functions: 



sin(4>) 

tan(<J>) 



* - -jU 3 

6 

* + 1* 3 



Therefore, if terms higher than third order are neglected: 



tan- 



: (4>)sin(4>) = f<t> * --g* 3 ) = + 2 f<l>-i+ 3 



- <t > 3 



and substituting into Equation (2.8) yields: 



GZ(4>) = GAT 




+ |rM<J> 3 = GM<j> + 



(-BM - -GM 

U 6 , 




Now for most ships, 3 BM>GM> which means that the leading cubic coefficient 
^ the GZ (<p) curve is positive and the curve concaves up as would a hardening 
spring. 

F. STATE-SPACE REPRESENTATION 

Combining Equations ( 2 . 2 ) and (2.4) and rearranging the terms, the equations 
of motion can be expressed in vector form: 
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Z = F( Z,6) 



(2.9) 



where z is the state vector [v r p <|>] T . 

Equation (2.9) is more conveniently expressed in matrix notation as is shown in 
Equations (2.10). The last matrix equation is derived from the kinematic roll rate. 



(mx c -NJ (I Z -NJ 0 -ty 




V 




/r 


(m-iy (mx G -Y f ) -mz G -Y^ 




r 




fs 


-mz G 0 (7,-tfp 




P 




fR 


0 0 0 1 




4> 




P. 



f Y = Nv+N^+N^vr 2 +(H r -mx G u)r+N rT / 3 +N vvr v 2 r+N^ + 

Y^+Y^+Y^+Y^rtf+Yw+Y, 

/* = V + ^wv v 3 + ^m v lr +(K r +m Z( M)r+K rTT r 3 + AT^vr 2 +^4) + 
^v 2 <}) +^v<}) 2 + ^r 2 <t> + *^r<t> 2 +K(6) +AT 0 -A GZ(<J>) 



(2.10a) 



(2.10b) 



(2.10c) 



(2.10d) 
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III. STEADY STATE REPRESENTATION 



A. LINEAR STEADY STATE SOLUTIONS 

With the above third order Taylor expansion, the equations of motion 
represented by Equations (2.10) are highly nonlinear, the solution of which depends 
upon the initial and equilibrium conditions. The steady state condition requires the 
acceleration terms ^ ^ p > 4 , to be equal to zero. Thus Equation (2. 10a) reduces to 



Equation (3.1) for steady state: 

fr 

fs 

h 

P. 



[ 0 ] 



(3.1) 



where f Y , f s , and f R are defined in Equations (2.10). 

The numerical solution to Equation (3.1) is an iterative process which depends 
upon a fairly accurate initial estimate of the root for convergence. This starting point 
can be determined by assuming an equilibrium state of v o =r o =p o =0 o =O and a very 

small initial rudder deflection ( 6 0 ^ 1 '). Application of these assumptions to Equation 

(3.1) yields the linear, steady-state, yaw-sway- roll equations of motion in Equation 

(3.2) . 
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(3.2) 



N v (N r -mx G u) 




V 




-N(6)-N 0 


Y v (Y r -mu) 




r 


= 


-Y(b)-Y 0 


K v (K r+ m Zc u) (K^-AGM) 




<f>. 




-K(b)-K 0 



Equation (3.2) is in the familiar matrix form Az = B, where z is the state vector 
with p=0. 

Various computer subroutines are commercially available which solve this 
simple system of equations for the state vector values, given the other parameters. 
The values for the state vector obtained in this manner are then used as the initial 
estimate in solving the nonlinear system of equations. 

B. NONLINEAR STEADY STATE SOLUTION 

Once the initial estimate for equilibrium is obtained, the system of nonlinear 
equations in Equation (3.1) can be solved numerically for the state vector nominal 
points for each rudder angle, 8. Here again, commercial computer subroutines are 
available. The Levenberg-Marquardt algorithm with analytic Jacobian was chosen as 
the solution method [Ref. 17]. Since p=0 in Equation (3.1), the system of 
nonlinear equations is reduced to three and the Jacobian was determined by 
Equations (3.3). 
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in 



(3.3a) 



Vy 


Vy 


Vy 


dv 


dr 


64> 


v s 


V s 


Vs 


dv 


dr 


d4> 


% 


v R 


V R 


dv 


dr 


d4>. 



Vy 

dv 






i. mt) 



dv 



^^N„vr*(N r -mx^)*3N m r 1 *N m v 2 *2N^*N^‘*^^ 

^‘2Y„vr*(Y r -mu)+3Y^Y^2Y^+Y„tf+^ 
V s „ .„ _x 1 _. a . 37(6) 



— = Y s + Y^v 2 +2Y V 4.v4>+Y r .r 2 +2Y rA .r 4> + 

64, ♦ ""♦ "W Y rri, "M> Y 04, 






dK(6) 

dv 



dr 



Y r ■“J'W + 3*> I+ J<wV^2*^r* + A-^ + 

It - ^v 2 4-2^^ + .2*^r4> <- - a(AGZ(<t> )) 



34> 



2j.p- ,,2 + 2 p - ,A±r a2x 3AT(6) 
64> 



34> 



(3.3b) 

(3.3c) 

(3.3d) 

(3.3e) 

(3.30 

(3-3g) 

(3.3h) 

(3.3i) 

(3-3j) 
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dN(6) m Q 



N(6) = F,(v,r) 






where Y(S) = F 2 (v,r) 
K(6) = F 3 (v,r) 



§jm m 



(3.3k) 



d4> 



mv m o 



a<j> 



and 



<3(AGZ(4>)) 

a<l> 



= A GMcos(4>)+fiMtan 2 (4>) 




(3.31) 



Since F„ F 2 , and F 3 are complicated functions of v and r (Equations (2.6) and 
(2.7)), the following approximation method was used for evaluation: 



Each term in the Jacobian is determined by evaluation at the previous nominal 



Section A above. 

The result of the numerical solution to the nonlinear equations (Equations 
(3.1)) is the steady state values for v, r, and <f> for each rudder angle under steady 
turning conditions and for a given GM and propeller speed. 

C. TYPICAL RESULTS 

Utilizing the solution method depicted in Sections A and B above, PROGRAM 
COUPLED listed in Appendix A, was developed to predict the steady state roll 



dF 

dx 




FQCp) - F(Q.99x 0 ) 
*0 ■ °-"*0 



(3.4) 



point with the initial nominal point determined from linear solutions as described in 
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angle, sway velocity and turning rate of a high speed container ship as a function of 
rudder angle. 

The design characteristics of a SR- 108 container ship are presented in Table 
IV. The hydrodynamic derivatives and coefficients used in the computer simulation 
are listed in Table V. 

For a base-line model, the hydrodynamic derivatives and coefficients were held 
constant in addition to the Froude number (hence the propeller revolutions), the 
transverse metacentric height, and speed loss ratio. Therefore, for the figures 
presented in this section: 



Figures 3.1 through 3.3 show the variation of sway velocity, yaw rate, and roll 
angle as a function of rudder angle during a steady turn. It can be seen that all 
steady state variables v, r and $ are highly nonlinear functions in 5 which is 
attributed to both the effect of the nonlinear terms in our equations of motion and 
the fact that the ship speed is reduced during the turn. Had a constant speed, linear 
model been utilized, v, r and (f> would appear as straight lines versus 8. 



Froude number 
Metacentric height 
Speed loss ratio 



GM = 0.3 m 
a = 0.6 



Fn = 0.3 
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TABLE IV. PRINCIPAL DIMENSIONS OF SR-108 CONTAINER SHIP [Ref. 9] 



Items 


Ship 


Model 


Hull Length B. P. 


L 


(m) 


175.00 


3.000000 


Breadth 


B 


(m) 


25.40 


0.435000 


Draft Fore 


Tp 


(m) 


8.00 


0.137100 


Aft 


t a 


(m) 


9.00 


0.154300 


Mean 


T 


(m 3 ) 


8.50 


0.145700 


Displacement Volume 




(m 3 ) 


21,222 


0.106860 


Height from keel to 










transverse metacenter 


KM 


(m) 


10.3900 


0.17810 


Height from keel to 










center of buoyancy 


BM 


(m) 


4.6154 


0.07912 


Block coefficient 


Cp 






0.55900 


Prismatic coefficient 


Cb 






0.58000 


Waterplane area 








0.68600 


coefficient 


Cw 






0.518 L 


Midship section 


On 






0.240 L 


coefficient 










LCB from forward 










Perpendicular 




Bilge keel 










Length 


L d 


(m) 


43.75 


0.7500 


Depth 




(cm) 


45.00 


0.7714 


Rudder Area 


Ar 


(m 2 ) 


33.0376 


0.009709 


Height 


H 


(m) 


7.7583 


0.133000 


Aspect Ratio 


A 






1.821900 


Area Ratio 


A R /L d 


(m) 




1/45.0 


Propeller 










Diameter 


D 


(m) 


6.533 


0.112 


Pitch Ratio 


P 






1.009 


Expanded Area Ratio 








0.670 


Boss Ratio 








0.180 


Number of blades 








5 
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TABLE V. HYDRODYNAMIC DERIVATIVES AND COEFFICIENTS [Ref. 9] 



A.) Hull Only 



Y„ 


-0.0116000 


N v 


-0.0038545 


K, 


0.00030260 


Y, 


0.0024200 


N r 


-0.0022200 


K 


-0.00006300 


Y, 


-0.0000630 


N, 


-0.0001424 


*4 


-0.00002100 


Y w 


-0.1090000 


N™ 


0.0014920 




0.00284300 


Y m 


-0.0405000 


N m 


0.0015600 


Kn 


0.00105650 


Y m 


0.0017700 


N m 


-0.0022900 


Kn 


-0.00004620 


Y w 


0.0214000 


N w 


-0.0424000 


K™ 


-0.00055800 


Y^ 


0.0460500 


N w , 


-0.0190580 


Kw, 


-0.00120120 


Y v ^ 


0.0030400 


N.* 


-0.0053766 


Km 


-0.00007930 


Y„, 


0.0093250 




-0.0038592 


K* 


-0.00024300 


Y t * 


-0.0013560 




0.0024195 


Km 


0.00003569 



B.) Propeller and Rudder 



a H 


0.237 


T 


1.090 


(l-w p ) 


0.816 


X H 


-0.480 


€ 


0.921 


Xr 


-0.500 


C RX 


0.710 


k 


0.631 


*P 


-0.526 


Zr 


0.033 


C Sr 


-0.156 


K t 


0.527-0.455J 


C pv 


0.000 


C 8rrr 


-0.275 


y 


0.088 v>0 
0.193 v<0 


C pr 


0.000 




1.960 


N P 


79.10 Fn=0.2 
118.64 Fn=0.3 
158.19 Fn=0.4 
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Figure 3.1: Sway Velocity Under Steady Turning 




RUDDER ANGLE (deg) 

Figure 3.2: Turning Rate Under Steady Turning 
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rv. STABILITY ANALYSIS 



A. COUPLED STEERING AND ROLL EQUATIONS OF MOTION 

Once the nominal points are determined for each rudder angle, Equation 
(2.10a) can be reduced to state-space representation (Equations (3.4)) for stability 
analysis. 



Bi = Az 



where 





V 




V 


i = 


r 




r 


P 


; z = 


P 






<t> 




4> 



B = 



(mx G -NJ (I Z -NJ 0 



On-YJ 


(mx 


a-*,) 


-mz G 


-n 


-mz G 




0 






0 




0 


0 




1 




Vy 


Vy 


Vy 


Vy 






dv 


dr 


dp 


34 > 






V s 


v 5 


Vs 


V s 




A = 


dv 


dr 


dp 


d <t> 






V R 


v R 


V* 


Vr 






dv 


dr 


dp 


dc t> 






dp 


dp 


dp 


dp 






dv 


dr 


dp 


64 > 





(4-la) 



(4-lb) 



(4.1c) 
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Matrix B can be thought of as the generalized mass matrix for the problem 
while matrix A is the Jacobian matrix of f(z,S) evaluated at each nominal point. For 
stability analysis purposes, the variables v, r, p, and <f> are understood to represent 
small deviations of the actual variables from their respective steady state values as 
computed in the previous chapter. The rationale behind Equations (4.1) is the 
introduction of small local perturbations in v, r, p, and 4> superimposed on their 
steady turning values. Lyapunov’s linearization theorem establishes that stability 
properties of a nominal point for a nonlinear system can, in general, be deduced 
from the stability properties of the corresponding linearized system. 

The terms in the matrix A are as defined in Equations (3.3) with the following 
additions: 



dfy 


Vs 


. V* . 0 


dp 


dp 


dp 


dp 


_ V 


= V = o 


dv 


dr 


84 > 




dp 


= 1 




dp 





Stability can be determined from the solution to the eigenvalue problem: 

XBz = Az 



(4.2) 



| A - = 0 
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If all eigenvalues, A, have negative real parts, the nominal point is asymptotically 
stable. If at least one eigenvalue is positive, the nominal point is unstable. The 
degree of stability of the nominal point is therefore related to the real part of the 
eigenvalues. The more negative they are, the faster the exponential convergence of 
solutions in time to the nominal point. The imaginary part of the eigenvalues 
characterizes the frequency content of the oscillatory behavior of the system 
response. 

A subroutine was included in the computer model to solve this eigenvalue 
problem as a function of rudder angle. 

B. DECOUPLED STEERING EQUATIONS 

In order to distinguish the characteristic contributions of the sway and yaw 
motions from those of roll, the steering equations (Equations (4.3)) were decoupled 
from the maneuvering equations (Equations (4.1)), resulting in Equations (4.3). 



\mx G -NJ 


VrW ' 




V 




fys 


0 n-YJ 


(mx G -Y.) 




r 




fss. 



frs = ( 4 - 3b ) 

/re ■ N y v*N ym v , *N m vr' 1 *(N r -mxji)r*N m r , *N mr v 2 r*N(. S) ( 4 - 3c ) 

This decoupling is based on the assumption that all cross coupling coefficients 
between roll and sway/yaw are zero, which is the usual approximation made for 
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surface ships where roll is studied independently from sway/yaw. One of the 
motivations for this study is to evaluate the degree of accuracy of such an 
approximation. 

The additional subscript, S, distinguishes the decoupled steering forces and 
moments from those of the fully coupled system. Following the same format as 
before, the state-space formulation for stability analysis is given in Equations (4.4). 






where 



B s = 






; z , j = 



(mx G -N,) (I t -N) 
(m-YJ 0 mx G -Y .) ) 



(4.4a) 



(4.4b) 






$fys 



YS 



dv dr 

$[ss dfss 

dv dr 



— = N +3N v 2 +N r 2 +2N vr+^^Q 

dv V vw vrr WT ^ 



$ 



=2N vrr vr+(N r -mx G u)+3N rrr r 2 +N ViT v 2 + 

Vs 



dN( 6) 



dr 

% - 2y„ r vrt(y,-m«)43y m r J *K„ T v 2 *®M 

or or 



(4.4c) 

(4.4d) 

(4.4e) 

(4.4f) 

(4-4g) 
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The approximation method of Equation (3.4) was used to determine rudder 
force and moment contributions in the Jacobian matrix A. 

The stability attributable to yaw and sway as a function of rudder angle is thus 
determined by the solution of the decoupled steering eigenvalue problem in 
Equations (4.5). 



^•s^s^s = 

or 



(4.5) 



C. DECOUPLED ROLL EQUATIONS 

By decoupling the steering equations, the roll equations become uncoupled 
producing the following relationship: 





-V 




p 




fflR 


0 


1 




, 4 >. 




.P . 



(4.6a) 



/« = + *(6) - A GZ(4>) 



(4.6b) 



The additional subscript, R, refers to the decoupled roll forces and moments 
and GZ (<f>) is as given in Equation (2.8). Equations (4.7) are the state-space stability 
analysis formulations. 
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where 



Z p — 



(4.7a) 



Br = 






• 


0 


1 


(4.7b) 


&RR 


^RR 




(4.7c) 


dp 


dp 






dp 


dp 






dp 


3<j> 







V, 



RR 



dp 



= 0 



&XR = ^ _ 3(AGZ(<j>)) 
0<t> ^ 3<J> 



(4.7d) 



Jp = i 
dp 



& = o 

3<1) 



The expression — — j s given in Equation (3.3). 

c5<t> 



(4.7e) 



Therefore, the stability attributable to roll as a function of rudder angle is 
attained by solution of the decoupled roll eigenvalue problem of Equations (4.8). 



^rBrZ# ArZz 



or 

I Ar-XrB r \ = 0 



(4.8) 
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D. TYPICAL RESULTS 



Continuing with the base-line model described in Chapter III Section C, the 
roots of Equations (4.2), (4.5) and (4.8) were obtained using PROGRAM 
COUPLED (Appendix A). Figure 4.1 is a plot of the real roots of the fully coupled 
maneuvering equations as a function of rudder angle. As can be seen, there are three 
distinct components which were classified as due to roll or due to steering (sway/yaw) 
by decoupling the steering and roll equations as discussed in Sections B and C above. 
Figures 4.2 through 4.4 show the three components of the fully coupled system 
matched with their respective decoupled solutions of Equations (4.5) for steering and 
Equations (4.8) for roll. 

In Figure 4.2, the decoupled roll real roots are constant for each rudder angle, 
whereas the coupled roots vary with an increasing stability to about five degrees of 
rudder. The upper steering component shown in Figure 4.3 indicates very little 
deviation in the coupled and decoupled roots, however, the lower steering 
component in Figure 4.4 shows a deviation in the two roots as the rudder angle is 
increased. The steering eigenvalues become increasingly more negative as the rudder 
angle is increased. This indicates that the steering system is dynamically more stable 
for a non-zero rudder angle, and as a result, a steering system that is designed to be 
stable for straight line motions will be even more stable for motions along curved 
reference paths. 

The imaginary roots for the decoupled equations are zero, indicating that all 
imaginary roots obtained in the solution of the coupled equations are due to roll 
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only. Figure 4.5 shows that for rudder angles less than 15°, the coupled model 
produces larger values than the decoupled model and for rudder angles greater than 
15°, the reverse is true. Comparing Figures 4.2 and 4.5, we can see that the imaginary 
part of roll response, or the roll frequency of oscillation, can be predicted fairly 
accurately from the decoupled roll equations. This is not the case for the real part 
where decoupling the equations results in a severe underestimation. 

Figure 4.6 is a root locus plot of the roll eigenvalues for both coupled and 
uncoupled models. As can be seen, the fully coupled model is more stable than the 
decoupled model. 

The roll damping caused by roll is shown in Figure 4.7. The decoupled roll 
analysis shows less damping for each rudder angle than does the fully coupled model. 
Furthermore, the damping increases significantly for large rudder angles. It can be 
seen that the actual roll damping is significantly larger than what is predicted from 
the decoupled model, which may underestimate the actual damping ratio by a factor 
of two to three. 
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Figure 4.1: Coupled Model Real Eigenvalues 
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Coupled Model and Decoupled Model Roll Eigenvalues 
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RUDDER ANGLE (deg) 

Figure 4.3: Coupled Model and Decoupled Model Upper Steering Eigenvalues 




RUDDER ANGLE (deg) 

Figure 4.4: Coupled Model and Decoupled Model Lower Steering Eigenvalues 
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RUDDER ANGLE (deg) 



Figure 4.5: Coupled Model and Decoupled Model Roll Eigenvalues 




Response 
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RUDDER ANGLE (deg) 

Figure 4.7: Coupled Model and Decoupled Model Damping Ratio 
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V. SENSITIVITY OF MATHEMATICAL MODEL 
TO VARIATION IN DESIGN PARAMETERS 



A. BACKGROUND 

Since actual ships rarely operate under the precise conditions assumed in the 
development of this mathematical model, changes in certain variables to measure the 
response sensitivity is prudent, as for any design study. The following sensitivity 
analysis has this aim in mind as well as to identify any inexplicable response as an 
area for further study. It must be remembered that certain assumptions were made 
at the outset which might need refinement as the design spiral continues. For 
instance, the roll damping term, K^, is assumed to be constant throughout the 
simulation, although experimental evidence [Ref. 18] suggests that roll damping 
coefficients are forward speed and motion amplitude dependent. By isolating a 
parameter, problem areas are enhanced so that a particular parameter can be the 
focus for analysis of the validity of the assumptions made or the modelling method 
used. 

If, on the other hand, one or more of the dependent variables do not show 
appreciable sensitivity, further simplifying assumptions could possibly be made when 
analyzing the more difficult problems of maneuvering in regular and irregular seas. 

In the following sections, sensitivity to speed loss during a turn, Froude number 
of approach spped, metacentric height, operating turn, and righting moment is 
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evaluated. For brevity, a limited number of plots are presented in each case. 
Appendix B contains additional plots for metacentric height, Froude number, and 
speed loss. 

B. SPEED LOSS DURING A TURN AS THE DESIGN PARAMETER 

The sensitivity of the model to changes in speed loss during the turn was 
investigated by holding transverse metacentric height (GM) and Froude number 
constant: 

GM = 0.3 m 
Fn = 0.3 

In determining the speed loss under steady turning conditions, it was assumed 
that the water depth was sufficient, i.e. >110 meters, to neglect this effect [Refs. 2, 
3]. However, it has been determined by several sources: Davidson (1944), Shiba 
(1960), Strom-Tejsen (1965), Eda and Crane (1965), among others [Ref. 2], that 
speed loss in steady turning is a function of the hull configuration (including rudder) 
and the turning diameter. Various other sources, [Refs. 3, 9, 10] consider rudder 
angle as an additional parameter in determining speed loss. In an effort to combine 
these effects, Figures 25 and 178 and Tables 29(a) and 29(b) from Principles of Naval 
Architecture [Ref. 2], Figure 11 of Eda and Falls [Ref. 3], and Figure 10 of Son and 
Nomoto [Ref. 9] were used to develop the following simple relationship to represent 
speed loss: 
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where 



— = 1 - a * 6 

V 

r o 

V = steady turning speed 
V 0 = approach speed 

a = speed loss ratio 
6 r = rudder angle in radians 
Variation of the speed loss parameter ranges from no speed loss (a=0) to maximum 

speed loss (a =0.73) during a turn with fixed controls. Results of these effects on 

sway velocity, yaw rate and roll angle may be viewed in Figures 5.1 through 5.3. The 

influence of speed loss on r is negligible, indicating that in all cases considered, the 

ship steady turning rate performance is not affected. On the contrary, v and 4> show 

a relatively significant modification due to speed loss. The sign reversal of the roll 

angle plot (Figure 5.3) for large rudder angles is rather uncommon, but it is, 

nevertheless, possible as a result of the nonlinearities that are present in the 

equations of motion for this model. 

The characteristic roots as a function of rudder angle were determined to 
establish stability trends. The real roots indicate three distinct components as 
previously discussed in Chapter IV, Section C. Both the upper and lower steering 
components demonstrate a significant increase in stability with rudder angle but little 
variance in stability due to speed loss. The roll component, however, indicates 
greater stability as speed loss decreases. We discovered earlier, that the non-zero 
values in the root locus were due solely to the roll contribution owing to the 
imaginary parts of the steering elements. The damping ratio of the roll contribution 
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RUDDER ANGLE (deg) 

Figure 5.1: Effect of Speed Loss on Sway Velocity 




Figure 5.2: Effect of Speed Loss on Yaw Rate 
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Figure 5.3: Effect of Speed Loss on Roll Angle 



imaginary parts of the steering elements. The damping ratio of the roll contribution 
is related to both real and imaginary parts of the eigenvalues. Figures 5.4 and 5.5, 
for the fully coupled and uncoupled models respectively, show this as a function of 
rudder angle. Both models exhibit qualitatively the same trend, with the uncoupled 
model showing a significant underestimation of the true damping ratio value. 

C. FROUDE NUMBER AS THE DESIGN PARAMETER 

In testing the model sensitivity to Froude number changes, the transverse 
metacentric height and speed loss ratio were kept constant: 

GM = 0.3 m 
a = 0.6 
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Figure 5.5: Effect of Speed Loss on Decoupled Model Damping Ratio 
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Froude numbers used to calculate the forward velocity prior to initiating the turn 
aswell as propeller revolutions were assumed to remain constant during the turn. For 
each Froude number and rudder angle, the sway velocity, yaw rate and roll angle 
were determined using the model program. Figures 5.6 through 5.8 at the end of this 
section depict graphic results. Note that v and r are not readily affected by change 
in Froude number, whereas <f> is greatly changed, increasing dramatically with Froude 
number. 

Using these values, the stability of the coupled response was determined. The 
real roots of the coupled maneuvering equations for the three Froude numbers are 
indistinguishable when plotted as in Figure 4.1, as are the roots for the lower 
steering components (Figure 4.2). The upper steering component exhibits nonlinear 
behavior through five degrees of rudder for all Froude numbers with increasing 
stability as both Froude number and rudder angle increase. The roll component was 
nonlinear throughout the range of rudder angles tried. Furthermore, as the Froude 
number increased, so did the range of rudder angles for which stability increased. 

The coupled model roll damping ratio remains relatively insensitive to Froude 
number (Figure 5.9), while the uncoupled model suggests a wide variation (Figure 
5.10). This demonstrates the incorrect conclusions which can be reached when using 
linear/uncoupled models that neglect the actual coupling of the lateral sway/yaw 
motions back into roll. 
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Figure 5.6: Effect of Fn on Sway Velocity 




Figure 5.7: Effect of Fn on Yaw Rate 
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Figure 5.9: Effect of Fn on Coupled Model Damping Ratio 
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D. TRANSVERSE METACENTRIC HEIGHT AS THE DESIGN PARAMETER 
For transverse metacentric height as the design parameter in the model, Froude 
number and speed loss ratio were held constant: 

Fn = 0.3 

a = 0.6 

Here again, the sway velocity, yaw rate and roll angle for each GM and rudder 
angle were computed using the model program. The results are represented in 
Figures 5.11 through 5.13. As indicated for the case of Froude number being the 
design parameter, v and r are affected very little by even an order of magnitude 
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change in GM. On the other hand, the increase in GM not only causes a large 
decrease in the magnitude of the roll angle, but in the range as well. 

In determining the influence of GM on the stability characteristics, it was found 
that the real roots plotted as a function of rudder angle showed no distinction 
between GM’s. When the components were separated out, only the upper steering 
and roll constituents exhibited significant variation. The upper steering response 
reflected a somewhat incongruous reaction, that is, as the GM increases, the stability 
due to steering decreases slightly. The roll component demonstrated a marked 
increase in stability with increase in GM. This result is evident in the root locus 
diagram of Figure 5.14. 

Conversely, the damping ratio did not change significantly, as shown in Figure 
5.15. Computation of the damping ratio from the decoupled model showed a much 
wider deviation among different GM’s. 

E. TRIM AS THE DESIGN PARAMETER 

For the investigation of a vessel in the trimmed condition, the transverse 
metacentric height, Froude number and speed loss ratio were held constant: 

GM = 0.3 m 
Fn = 0.3 
a = 0.6 
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Figure 5.12: Effect of GM on Yaw Rate 
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Figure 5.14: Effect of GM on Coupled Model Root Locus 
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Figure 5.16: Effect of GM on Decoupled Model Damping Ratio 
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The base-line model chosen for this study, and therefore the hydrodynamic 
derivatives and coefficients derived from model tests, were for a design trim of 0.5 
meters by the stern and a mean draft of 8.5 meters for the actual ship (refer to Table 
IV for model parameters). It has been shown [Ref. 11] through linear theory that 
these derivatives and coefficients change for the trimmed condition by the following 
semi-empirical relations: [Ref. 15] 

W = 

N v (t) = 

w = 

N r ( T) = 



where 

x = 

d n = 

T v (0) = 

The application of these relationships are somewhat questionable, but in the 
absence of model tests, they will suffice to indicate trends for further study. 

A further modification in the model was made to the z R parameter (Table III) 
or the z coordinate of the point on which the rudder force acts: 



Y v (0) 
N v (0) 
Y r (0) 
N r ( 0) 



1 + 



1 - 



2x 
3 d 

m 

0.27 x 



l d 

v m 



l + 



l + 



0.8 x 



m 

0.3 x 



/ - Ml 

v T v (0) 



+ for trim by stem 
- for trim by bow 
design mean draft 
Y at design condition 
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Figures 5.17 through 5.19 depict the effects of — = ±0.2 which equates to a 
significant change in operating condition. The sensitivity of v and r to the previous 

parameters was practically negligible with the exception of the speed loss ratio. This 
is not the case for trim where turning rate is improved for trim by the bow and 
decreased for trim by the stern. It should be noted that this is contrary to linear 
theory results. For sway velocity at small angles of rudder the variation is small, but 
for large rudder angles the trim by stern curve deviates considerably from the other 
two conditions. The difference is more pronounced in the roll angle curves (Figure 
5.19) where trim by the bow generates higher angles during the turn. On the other 
hand, trim by the stern results in smaller roll angles with a roll angle sign reversal 
beyond a certain rudder angle. This is due to the fact that for large rudder angles, 
the hydrodynamic moment exerted by the rudder force dominates over the fluid 
forces on the hull. 

The root locus diagram (Figure 5.20) for the roll stability analysis is revealing 
in that very different trends can be observed. Trim by the stern is considerably more 
stable than by the bow for the entire range of heel angle, although at larger rudder 
angles, the stability is less than at the design condition. This is possibly qualitatively 
explained by noticing in Figures 5.21 and 5.22, that the damping ratio for trim by the 
stern increases sharply precisely at the inflection point of the roll angle as shown in 
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Figure 5.19. Here again, it must be pointed out that the damping ratio’s variance 
between the coupled model and the decoupled model is at least a factor of two. It 
can be seen from both coupled and decoupled models that at the inflection point of 
the roll angle, the damping ratio reaches its peak value. This maximum value of the 
damping ratio exists for the other trim conditions as well, but its actual location is 
shifted outside the range of rudder angle variations considered here. 

F. RIGHTING ARM CURVE AS THE DESIGN PARAMETER 

The sensitivity of the model to changes in the righting arm as a function of the 
roll angle was tested by holding the transverse metacentric height, Froude number 
and speed loss constant: 



GM = 1.0 m 
Fn = 0.3 
a = 0.6 



In general, the righting arm curve, GZ (4>), is taken to be a very simplified 
version of Equation (2.8), namely GZ(<£)=GM*<£. This linear righting arm gives 
good approximations for small angles of heel and is derived much the same as was 
done for Equation (2.8) in Chapter II, Section E. For the purpose of this 
investigation, the simplified GZ(</>) was left in its nonlinear form, 
GZ(<£)=GM*sin(<£). Figure 5.23 illustrates the two GZ (<f>) curves used. Note that 
GZ1 (representing Equation (2.8)) develops a much larger righting arm at increased 
angles of heel than does GZ2 (representing GM*sin(<£)). 
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Figure 5.17: Effect of Trim on Sway Velocity 




RUDDER ANGLE (deg) 

Figure 5.18: Effect of Trim on Yaw Rate 
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Figure 5.19: Effect of Trim on Roll Angle 




REAL PART 

Figure 5.20: Effect of Trim on Coupled Model Root Locus 
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RUDDER ANGLE (deg) 

Figure 5.22: Effect of Trim on Decoupled Model Damping Ratio 
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The effects of these different righting arm curves on the model for v, r, and <£, 
are shown in Figures 5.24 through 5.26. Sway velocity and yaw rate appear to be 
unaffected by this parameter change in contrast to roll angle, which shows relatively 
high sensitivity, as expected. 

Figures 5.27 (coupled model) and 5.28 (decoupled model) indicate the same 
trend in damping ratio sensitivity as the previous parameters, specifically, the 
decoupled model underestimates c by a factor of two. 

G. SUMMARY OF RESULTS 

A qualitative summary of results is presented in Table VI below. All 
comparisons are to the base-line configuration of: 1 

GM = 0.3 m 
Fn = 0.3 
a = 0.6 



1 For GZ sensitivity GM = 1.0 m. 
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Figure 5.23: The Two Righting Arm Curves 




Figure 5.24: Effect of GZ (<t>) on Sway Velocity 
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Figure 5.25: Effect of GZ (<£) on Yaw Rate 




RUDDER ANGLE (deg) 

Figure 5.26: Effect of GZ (<f>) on Roll Angle 



62 





RUDDER ANGLE (deg) 



Figure 5.28: Effect of GZ($) on Decoupled Model Damping Ratio 
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Notation in the table refers to the largest noted deviations from the model 



according to the following subjective scale: 



Sensitivity 


Symbol 


Standard 


Negligible 


N 


< i% 


Low 


L 


1% - 5% 


Moderate 


M 


5% - 20% 


High 


H 


> 20% 


Mixed 


X 


=> large change 
over the range 



Additionally, A, B and C refer to rudder angle range for the noted deviations such 
that: 



0’ < A z 10° 

10 ° < B z 20° 

20° < C z 30° 
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TABLE VI. QUALITATIVE SENSITIVITY SUMMARY 



Parameter 


8 


a 


Fn 


GM 


Trim 


GZ{4>) ”1 


V 


A 


M 


N 


N 


L 


N 




B 


M 


N 


N 


M 


N 




C 


H 


N 


N 


H 


N 


r 


A 


N 


L 


N 


L 


N 




B 


N 


L 


N 


L 


N 




C 


L 


L 


N 


L 


N 


4> 


A 


M 


H 


H 


X 


H 




B 


M 


H 


H 


H 


H 




C 


H 


H 


M 


H 


M 


coupled 


A 


M 


M 


X 


X 


X 


damping 


B 


L 


M 


H 


X 


H 




C 


X 


X 


X 


X 


H 


uncoupled 


A 


L 


H 


X 


X 


X 


damping 


B 


M 


H 


H 


X 


H 




C 


H 


H 


H 


X 


H 



This reference table shows clearly that the sway velocity and turning rate are 
relatively insensitive to Froude number, metacentric height and righting moment, 
while being quite sensitive to the speed loss and trim parameters. Roll and damping 
ratio show the greatest response for all the parameters considered. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

This thesis developed a mathematical model for the coupled sway, yaw and roll 
equations of motion with rudder deflection as the independent variable. Additionally, 
various components of the hydrodynamic forces and moments, e.g., rudder to hull 
interaction, were modelled using existing research formulation. Deviation of this 
study from current solution methods was in the manner of evaluating the Taylor 
expansion for the steady state Jacobian matrix. The equilibrium position for each 
rudder angle was found to deviate significantly from the commonly used zero state, 
even for rudder angles of less than 10°. In fact, the roll angle and sway velocity 
nominal points are nearly maximum at 10° of rudder. 

Comparison of stability results for the coupled lateral motion model to those 
of models representing the steering equations of motion decoupled from the roll 
equations of motion showed: 



1. The directional stability characteristics signified by the real part of the 
complex eigenvalues are not affected significantly by the coupling of roll into 
sway and yaw. This indicates that sway and yaw can be treated separately 
from roll when studying directional stability. 

2. Roll stability represented by the magnitude of the damping ratio is 
considerably affected by the coupling of sway and yaw into roll. The fact that 
the damping ratio is underestimated by a factor of two by the uncoupled 
treatment of sway and yaw, indicates ship design by this means is a significant 
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over-design that bears more attention, especially with the economic 
constraints of today’s Navy. 

3. The natural frequency is not significantly affected by coupling. This indicates 
that while maneuvering in a seaway studies require the coupling of the lateral 
motions, the induced roll from the static righting moment can be studied from 
decoupled systems. 

4. Steering stability increases with increasing rudder deflection. This may 
indicate that if a ship is directionally unstable for straight line motion, it 
becomes directionally stable during a turn through the coupling of the lateral 
motions. 

The sensitivity of certain design parameters to the coupled system was 
investigated with the following determinations: 



1. The parameters tested were speed loss during a turn, transverse metacentric 
height, Froude number of initial forward speed, trim and righting moment. Of 
these parameters, trim proved to be the most affected. This indicates that 
during all phases of design, operating trim must be considered in addition to 
the design trim. 

2. Sway velocity and turning rate were relatively insensitive for all parameters 
except for trim and speed loss. 

3. Roll angle and damping ratio of all the parameters were significantly affected. 

4. The coupled model damping ratio showed significantly lower sensitivity than 
the decoupled model for small rudder angles in all parameters. 



B. RECOMMENDATIONS 

The results produced in this thesis were as a consequence of experimentally 
reported nonlinear hydrodynamic coefficients and design parameters of a SR- 108 
container ship [Ref. 10]. Additional studies using similar data are required. 
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Utilizing the method of determining equilibrium points for each rudder 
position, design a control law for an autopilot and determine the effects of the 
coupled lateral motions on antopilot dynamics. 

The sharp change in roll angle ( <f > ) for small values of rudder angle (5) suggests 
that this ship would be an excellent candidate for rudder roll stabilizing techniques. 

Finally, as indicated in Chapter I, Section B, use the coupled sway, yaw and roll 
equations of motion and build a mathematical model which treats wave effects in 
regular and irregular seas. Once this model is developed, analyze maneuvering in 
waves from the standpoint of biased roll oscillations. Biased roll oscillations are 
meant as roll about non-zero heel angle. 
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APPENDIX A 



*********************************************************************** 
★ * 

★ 
★ 



PROGRAM COUPLED 

This program solves the coupled yaw-sway- roll ship maneuvering 
equations. In addition, the yaw-sway steering equations are de- 
coupled from the roll equations and solved in order to compare 
results with the fully coupled results. The program is broken 
into several subroutines, both original and standard commercial 
programs . 



★a********************************************************************* 



★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★a************************ 



★ ★ 

* The INCLUDE statement is used to pass all the variables between * 

* the many subroutines. * 

* ★ 



INCLUDE ’ HYDCOMN.for’ 

COMMON /WORKSP/RWKSP 
REAL RWKSP(5000) 

REAL*8 VDEL (300) ,V(300) ,R(300) ,RPHI(300) 

EXTERNAL FCN , LSJAC , DNEQNJ 

★ 

* Open the input and output files 

★ 

OPEN ( 10 , FILE =’C0EF. INP’ , status= ’ old ’ ) 

0PEN(15, FILE =’VAR.0UT’ ,status=’new’ ) 

0PEN(20, FILE =’ NLIN.DAT’ , status= ’ new’ ) 

0PEN(25 , FILE = ’ EIGVAL1 . DAT’ , status= ’ new’ ) 

0PEN(30 , FILE = ’ EIGVAL2 . DAT’ , status= ’ new ’ ) 

OPEN (26, FILE = ’ EIGVALL . DAT ’ , status= ’ new’ ) 

0PEN(27, FILE =’ EIGVALp.DAT’ ,status=’new’ ) 

0PEN(35 , FILE =’ ZETA.DAT’ , status= ’ new’ ) 

OPEN (36 , FILE =’ OMEGA . DAT’ , status= ’ new’ ) 

0PEN(37, FILE = ’ZROLL .DAT’ , status= ’ new’ ) 

0PEN(38, FILE = ’ RUDDER . DAT ’, status= ’ new’ ) 

PI=4 . 0D0*DATAN ( 1 .0D0) 

ERRREL=0 . 0001 
I TM AX =200 

★ 

This subroutine reads the hydrodynamic derivative and coefficient 

* values from COEF.INP and writes pertinent values in VAR. OUT. 

CALL INPUT 
DEL=DEL* (PI / 1 80 . 0D0) 

* 

* This subroutine sets up the linear EOM’s in matrix form: 

* A x = B 

* 

CALL LINEAR 

* 

* This is an IMSL MATH/LIBRARY Subroutine [Ref. 19} 
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Purpose: Solve a real general system of linear equations with 

iterative refinement 



Arguments : 

N -Number of equations. (Input) 

A -N by N matrix containing the coefficients of the linear 
system. (Input) 

LDA -Leading dimension of A. (Input) 

B -Vector of length N containing the right-hand side of the 
linear system. (Input) 

IPATH -Path indicator, if = 1 the system A*x=B is solved 
X -Vector of length N containing the solution to the linear 
system. (Output) 

CALL DLSARG ( N , A , 3 , B , 1 , XGUESS ) 

WRITE(15,*) ’THE SOLUTION TO THE LINEAR SYSTEM OF EQUATIONS IS’ 
WRITE (15,*)’ ’ 

WRITE ( 15 , * ) ’V0=’ ,XGUESS(1) 

WR ITE ( 1 5 , * ) ’ R0= ’ , XGUESS ( 2 ) 

WRITE(15,*) ’PHIO=’ ,XGUESS(3) 

WRITE (15,*)’ ’ 

PH I0=XGUESS ( 3 ) * ( 1 80 . 0D0 / P I ) 

WRITE (15,*)’ PHIO IN DEGREES =’,PHI0 
WRITE ( 15, *) ’ ’ 

Set the rudder limits for steady turning 

DELMIN=0.01D0 
DELMAX=30 . 0D0 

This Do-Loop calculates the values of sway velocity, yaw rate and 
roll angle for each rudder angle. 

DO 100 IDEL=1 , IMAX 
WRITE (*,*) IDEL , IMAX 

RDEL=DELMIN+(DELMAX-DELMIN) *REAL( ( IDEL-1 ) / (IMAX-1 ) ) 

DEL=RDEL* ( PI / 1 80 .0D0 ) 

UND=1 .ODO-SLR* (DABS (DEL) ) 

This is an IMSL MATH/LIBRARY Subroutine [Ref. 17] 

Purpose: Solve a system of nonlinear equations using the 

Levenberg-Marquardt algorithm with a user-supplied 
Jacobian . 

Arguments : 

FCN - User- supplied SUBROUTINE to evaluate the system of 
equations to be solved. The usage is 
CALL FCN (X, F, N), where 

X - The point at which the functions are 

evaluated. (Input) 

X should not be changed by FCN. 

F - The computed function values at the point X. 

(Output) 

N - Length of X, F. (Input) 

FCN must be declared EXTERNAL in the calling program. 
LSJAC - User- supplied SUBROUTINE to evaluate the Jacobian at a 
point X. The usage is 
CALL LSJAC (N, X, FJAC) , where 



* 

★ 

★ 

★ 

* 

★ 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

★ 

* 



★ 



* 

★ 

★ 



N - Length of X. (Input) 

X - The point at which the function is evaluated. 

( Input) 

X should not be changed by LSJAC. 

FJAC - The computed N by N Jacobian at the point X. 
(Output) 

LSJAC must be declared EXTERNAL in the calling program. 
ERRREL - Stopping criterion. (Input) 

The root is accepted if the relative error between two 
successive approximations to this root is less than 
ERRREL. 

N - The number of equations to be solved and the number 

of unknowns. (Input) 

ITMAX - The maximum allowable number of iterations. (Input) 

The maximum number of calls to FCN is ITMAX* (N+1). 
Suggested value = 200. 

XGUESS - A vector of length N. (Input) 

XGUESS contains the initial estimate of the root. 

X - A vector of length N. (Output) 

X contains the best estimate of the root found by 
DNEQNJ . 

FNORM - A scalar which has the following value, 

F(1 )**2+. . .+F(N)**2 at the point X. (Output) 



CALL DNEQNJ ( FCN , LSJAC , ERRREL , N , ITMAX , XGUESS , X , FNORM ) 

PH I =X ( 3 ) * 1 80 . 0D0 / P I 
WRITE(20,50)RDEL,X(1) ,X(2) ,PHI 
50 FORMAT (4E15 . 5) 

Calculate the rudder forces and moments: Equations 2.6 and 2.7 



UP=UND* (WP+TAU* ( (X ( 1 ) +XP*X (2) ) **2+CPV*X( 1 ) +CPR*X(2) ) ) 
JP=(U*UP) / (NP*D) 

KT=KTA+ (KTB*JP) 

UR=UP*EPSILON* ( (1 .0D0+8.0D0*K*KT) / (PI*(JP**2) ) )**0.5 
VR=GAMMA*X( 1 ) +CDR*X (2) +CDRRR* (X (2) **3) +CDRRV*X (1)*(X(2)**2) 
VD=(UR**2+VR**2) **0.5 
ALPHAR=DEL+DATAN(VR/UR ) 

FLAM= (6 . 13D0* LAMBDA) / (LAMBDA+2 . 25D0) 

FNN=FLAM* (AR/(L**2))*((VD)**2) *DSIN (ALPHAR ) 



YD=- (1 .ODO+AH) *FNN*DCOS(DEL) 
ND=- (XR+AH*XH) *FNN*DCOS(DEL) 
KD=(1 .ODO+AH) *ZR*FNN*DCOS(DEL) 
WRITE (38, *)ND,YD,KD 

Rudder angle in degrees 



VDEL(IDEL)=DEL 
V ( IDEL) =X ( 1 ) 

R ( I DEL ) =X ( 2 ) 



Rudder angle in radians 

RPHI ( IDEL) =X (3) 

XGUESS (1 )=X( 1 ) 
XGUESS(2) =X (2) 
XGUESS(3) =X (3) 

100 CONTINUE 
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This subroutine solves the eigenvalue problem |A - lambda*B|= 0 
for the fully coupled yaw-sway-roll equations, using the values 
for v, r, and phi obtained above. 

CALL EIGTH1 (VDEL,V,R,RPHI) 

This subroutine solves the eigenvalue problem |As - lambda*Bs|= 0 
for the steering yaw-sway equations. 

CALL EIGTH2(VDEL,V,R) 

This subroutine solves the eigenvalue problem |Ar - lambda*Br|= 0 
for the roll equations. 

CALL EIGTH3(VDEL,Rphi) 

STOP 

END 



SUBROUTINE INPUT 
INCLUDE ’HYDCOMN.for’ 

The numbers in the variables such as N122 stand for: 



READ(10,*)N0,N1 ,N2 , N3 , N1 1 1 , N1 22 ,N222 , N1 1 2 ,N1 13 , N133 ,N223 



& D , CDR , CDRRR , CDRRV , GAMMA 
Calculate the nondimensional mass 
M=W/ (0 .5D0* (L**3) ) 

WRITE(*,*) ’ENTER THE GM AND THE SPEED LOSS RATIO’ 

READ (*,*)GM, SLR 
WRITE (15,*) ’GM=’ ,GM 
WRITE (15,*)’ ’ 

WRITE (15,*) ’THE SPEED LOSS RATIO IS=’,SLR 

Adjusts the z coordinate of CG based on the GM 

ZGM=DABS ( 0 . 3D0 - GM ) / L 
ZG=0. 00078 -ZGM 

WRITE(15,*) ’THE SHIP DISPLACEMENT (IN CUBIC METERS) IS W=’,W 
WRITE ( 15 , * ) ’ ’ 

WRITE(15,*) ’THE SHIP NONDIMENSIONAL MASS IS M=’,M 
WRITE (15,*) ’ ’ 

XG=XG/L 

WRITE(15,*) ’L=’ ,L 
WRITE (15,*) ’XG=’ , XG 
WR ITE ( 1 5 , * ) ’ ZG= ’ , ZG 
WRITE(15,*) ’ ’ 

WRITE (*,*) ’ENTER THE INITIAL RUDDER ANGLE IN DEGREES’ 

READ( * , * )DEL 

WRITE(15,*) ’THE INITIAL RUDDER ANGLE DEL=’,DEL 
WRITE(*,*) ’ENTER THE FROUDE NUMBER AND THE IMAX’ 

READ ( * , * ) FN , IMAX 



1 implies 

2 implies 

3 implies 

4 implies 



v 

r 



phi 

delta 




WRITE (15,*) ’THE FROUDE NUMBER IS FN= ’ , FN 
Calculate the initial forward speed 
U=FN* ( (9 . 8D0*L ) **0 . 5D0) 

WRITE(15,*) ’THE FORWARD VELOCITY IS U=’,U 
WRITE (15,*) ’ ’ 

WRITE(*,*) ’ENTER THE PROPELLER SPEED IN RPMS FOR THE FN’ 
READ ( * , * ) NP 

WRITE(15,*) ’THE PROPELLER SPEED IN RPMS IS NP=’,NP 
Change rpm’s to radians per sec 
NP=NP*(PI/30.0D0) 

The weight component is nondimensionalized such that WG*GM is 
dimensionless 

WG= (W*9 . 81 DO) / (0 . 5D0* ( L**3) * (U**2) ) 

WRITE(15,*) ’THE NONDIM COMPONENT OF DISPLACEMENT IS WG=’,WG 

RETURN 

END 

SUBROUTINE LINEAR 

INCLUDE ’ HYDCOMN .f or ’ 

INTEGER IROW, ICOL 
DO 100 IR0W=1 ,N 
B( IROW) =0 .ODO 
DO 90 IC0L=1 ,N 
A( IROW, ICOL) =0 .ODO 
90 CONTINUE 
100 CONTINUE 

Equation 3.2 

A ( 1 , 1 ) =N1 

A ( 1 ,2)=N2- (M*XG*UND) 

A( 1 ,3)=N3 
A(2,1 )=Y1 
A(2 ,2) =Y2- (M*UND) 

A(2 ,3) =Y3 
A(3 , 1 ) =K1 

A(3 , 2) =K2+ (M*ZG*UND) 

A(3,3)=K3- (WG*GM) 

FNL=( (6. 13D0*LAMBDA) / (LAMBDA+2 . 25D0) )*(AR/L**2) 

N4=- (XR+ (AH*XH) ) *FNL 
Y4=- ( 1 .ODO+AH) *FNL 
K4=( 1 .ODO+AH) *ZR*FNL 

B(1 )=- (N4*DEL) -NO 
B(2) = - (Y4*DEL) -YO 
B(3)=- (K4*DEL) -KO 
RETURN 
END 

SUBROUTINE FCN(Z,Q,NDIM) 

INCLUDE ’HYDCOMN. for’ 

INTEGER NDIM 

DIMENSION Z(NDIM) ,Q(NDIM) 
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Equation 2.7 



★ 



UP=UND* (WP+TAU* ( (Z(1 )+XP*Z(2) ) **2+CPV*Z(1 )+CPR*Z(2) ) ) 

JP= (U*UP) / (NP*D) 

KT=KTA+ (KTB*JP) 

UR=UP*EPSILON*( (1 .0D0+8.0D0*K*KT)/(PI*(JP**2) ) )**0.5 
VR=GAMMA*Z( 1 ) +CDR*Z(2) +CDRRR* (Z(2) **3) +CDRRV*Z( 1 ) * (Z(2) **2) 
VD=(UR**2+VR**2) **0.5 
ALPHAR=DEL+DATAN(VR/UR) 

FLAM= (6.1 3D0*LAMBDA) / ( LAMBDA+2 . 25D0 ) 

FNN=FLAM* (AR/ (L**2) ) * ( (VD) **2) *DSIN(ALPHAR) 

★ 

* Equation 2.6 

★ 

Y4=- (1 .0D0+AH) *FNN*DCOS(DEL) 

N4=- (XR+AH*XH) *FNN*DC0S (DEL) 

K4= ( 1 . ODO+AH) *ZR*FNN*DCOS(DEL) 

* 

* Equation 2.10b 

* 

C1=N1*Z(1)+N111*(Z(1)**3)+N122*Z(1)*(Z(2)**2) 

C2=N222* (Z(2)**3)+N112*(Z(1)**2)*Z(2)+N3*Z(3) 

C3=N113*(Z(1)**2)*Z(3)+N133*Z(1)*(Z(3)**2)+(N2- (M*XG*UND) )*Z(2) 
C4=N233*Z(2) * (Z(3) **2) +N4+N0+N223* (Z(2) **2) *Z(3) 

★ 

Q ( 1 ) =C1 +C2+C3+C4 

★ 

* Equation 2.10c 

* 

D1=Y1*Z(1)+Y111*(Z(1)**3)+Y122*Z(1)*(Z(2)**2)+(Y2- (M*UND) ) *Z(2) 
D2=Y222* (Z(2)**3)+Y112*(Z(1)**2) *Z (2)+Y3*Z(3) 

D3=Y1 13*(Z(1)**2)*Z(3)+Y133*Z(1)*(Z(3) **2) +Y223* (Z(2)**2)*Z(3) 
D4=Y233*Z ( 2 ) * ( Z ( 3 ) * * 2 ) +Y4+Y0 

* 

Q ( 2 ) =D1 +D2+D3+D4 

★ 

* Equation 2.9 

★ 

c GZ=GM*DSIN (Z(3) ) +0 .5*BM* (DTAN (Z(3) ) ) **2*DSIN(Z(3) ) 

GZ=GM*DSIN(Z(3) ) 

* 

* Equation 2.10d 

★ 

El =K1 *Z(1)+K111*(Z(1) **3) +K1 22*Z( 1 ) * (Z(2) **2) + (K2+M*ZG*UND) *Z(2) 
E2=K222* (Z(2) **3) +K1 12*(Z(1)**2)*Z(2) +K3*Z(3) 

E3=K1 13*(Z(1)**2) *Z(3) +K133*Z( 1 ) * (Z(3) **2) +K223* (Z(2)**2)*Z(3) 
E4=K233*Z (2) * (Z(3) **2)+K4+K0- (WG*GZ) 

★ 

Q(3)=E1 +E2+E3+E4 

★ 

RETURN 

END 

* 

SUBROUTINE LSJAC(NDIM ,Z,QJAC) 

INCLUDE ’HYDCOMN.for’ 

INTEGER NDIM 

DIMENSION Z(NDIM) , QJAC ( NDIM , NDIM ) 

★ 

* Equation 2.7 
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UP=UND* (WP+TAU* ( (Z(1)+XP*Z(2) ) **2+CPV*Z ( 1 ) +CPR*Z(2) ) ) 

V9=Z ( 1 ) *0 .99D0 
R9=Z ( 2 ) *0 . 99 DO 

UPV9=UND* (WP+TAU* ( (V9+XP*Z(2) ) **2+CPV*V9+CPR*Z(2) ) ) 

UPR9=UND* (WP+TAU* ( (Z( 1 ) +XP*R9) **2+CPV*Z ( 1 ) +CPR*R9) ) 
JP=(U*UP)/(NP*D) 

JPV9= (U*UPV9) / (NP*D) 

JPR9=(U*UPR9)/(NP*D) 

KT=KTA+ (KTB*JP ) 

KTV9=KTA+ ( KTB* JPV9 ) 

KTR9=KTA+ (KTB*JPR9) 

UR=UP*EPSILON* ( ( 1 .0D0+8 .ODO*K*KT) /(PI*(JP**2)))**0.5 
URV9=UPV9*EPSIL0N* ( ( 1 .0D0+8 .0D0*K*KTV9) / ( PI* ( JPV9**2 ) ) ) **0 . 5 
URR9=UPR9*EPSIL0N* ( (1 .0D0+8.0D0*K*KTR9) / ( PI* ( JPR9**2) ) ) **0 .5 
VR=GAMMA*Z ( 1 ) +CDR*Z(2) +CDRRR* (Z(2 ) **3) +CDRRV*Z( 1 ) * (Z (2 ) **2) 
VRV9=GAMMA*V9+CDR*Z(2) +CDRRR* (Z( 2 ) **3) +CDRRV*V9* (Z ( 2 ) **2) 
VRR9=GAMMA*Z ( 1 ) +CDR*R9+CDRRR* (R9**3)+CDRRV*Z( 1 ) * (R9**2) 

VD= (UR**2+VR**2) **0 .5 
VDV9= (URV9**2+VRV9**2) **0.5 
VDR9= (URR9**2+VRR9**2) **0.5 
ALPHAR=DEL+DATAN ( VR / UR ) 

ALPHARV=DEL+DATAN(VRV9/URV9) 

ALPHARR=DEL+DATAN(VRR9/URR9) 

FLAM= (6.1 3D0* LAMBDA ) / ( LAMBDA+2 . 25D0 ) 

FNN=FLAM* ( AR/ (L**2))*((VD)**2) *DSIN ( ALPHAR) 

FNNV9=FLAM* (AR/ (L**2) ) * ( (VDV9) **2) *DSIN (ALPHARV) 

FNNR9=FLAM* (AR/ (L**2) ) * ( (VDR9) **2) *DSIN(ALPHARR) 

★ 

Y4=* ( 1 .ODO+AH) *FNN*DCOS(DEL) 

Y4V9=- (1 .0D0+AH)*FNNV9*DC0S(DEL) 

Y4R9=- (1 .0D0+AH)*FNNR9*DC0S(DEL) 

★ 

* Equation 3.4 

★ 

PY4V=(Y4-Y4V9) / (Z ( 1 ) -V9) 

PY4R=(Y4-Y4R9) / (Z(2) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS(DEL) 

N4V9=- (XR+AH*XH)*FNNV9*DC0S(DEL) 

N4R9=- (XR+AH*XH) *FNNR9*DC0S(DEL) 

PN4V= ( N4 - N4V9 ) / ( Z ( 1 ) - V9 ) 

PY4R=(Y4-Y4R9) / (Z(2) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS(DEL) 

★ 

* Equation 3.4 

★ 

N4V9=- (XR+AH*XH) *FNNV9*DC0S(DEL) 

N4R9=- (XR+AH*XH) *FNNR9*DC0S(DEL) 

PN4V= ( N4 - N4V9 ) / ( Z ( 1 ) - V9 ) 

PN4R=(N4-N4R9) / (Z(2) -R9) 

K4=(1 .ODO+AH) *ZR*FNN*DCOS(DEL) 

K4V9= ( 1 .ODO+AH ) *ZR*FNNV9*DC0S ( DEL ) 

K4R9= ( 1 .ODO+AH) *ZR*FNNR9*DC0S(DEL) 

★ 

* Equation 3.4 

★ 

PK4V= (K4-K4V9) / ( Z ( 1 ) - V9 ) 

PK4R=(K4-K4R9) / (Z(2) -R9) 

★ 

* Equation 3.3b 

★ 

C11=N1+3*N111*(Z(1) **2) +N122* (Z (2) **2) +2 .0D0*N1 12*Z ( 1 ) *Z (2) 
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* 



c 

c 

c 



* 



C12=2.0D0*N113*Z(1)*Z(3)+N133*(Z(3)**2) 

QJAC (1,1) =C1 1 +C1 2+PN4V 
Equation 3.3c 

C21 =2 .0D0*N1 22*Z( 1 ) *Z(2) + (N2-M*XG*UND) +3 .0D0*N222* (Z (2) **2) 
C22=N1 1 2* (Z ( 1 ) **2)+2 .0D0*N223*Z (2) *Z (3) +N233* (Z(3) **2) 

QJAC ( 1 , 2 ) =C21 +C22+PN4R 

Equation 3.3d 

C31=N3+N1 13* (Z(1)**2)+2.0D0*N133*Z(1) *Z (3) +N223* (Z (2) **2) 
C32=2.0D0*N233*Z(2) *Z(3) 

QJAC(1 ,3)=C31+C32 

Equation 3.3e 

D11=Y1+3.0D0*Y111*(Z(1 )**2)+Y122*(Z(2)**2)+2.0D0*Y112*Z(1)*Z(2) 
D12=2.0D0*Y113*Z(1)*Z(3)+Y133*(Z(3)**2) 

QJAC (2,1 )=D1 1 +D12+PY4V 

Equation 3.3f 

D21=2.0D0*Y122*Z(1)*Z(2)+(Y2-M*UND)+3.0D0*Y222*(Z(2)**2) 

D22=2 .0D0*Y223*Z (2) *Z (3) +Y233* (Z(3)**2)+Y112*(Z(1)**2) 

QJAC (2,2) =D21 +D22+PY4R 

Equation 3.3g 

D31 =Y3+Y1 13* (Z ( 1 ) **2) +2 .0D0*Y133*Z ( 1 ) *Z(3) +Y223* (Z(2) **2) 
D32=2.0D0*Y233*Z(2)*Z(3) 

QJ AC ( 2 , 3 ) =D31 +D32 

Equation 3.3h 

El 1 =K1 +3 ,0D0*K1 1 1 * (Z ( 1 ) **2) +K1 22* (Z(2) **2) +2 .0D0*K1 1 2*Z ( 1 ) *Z(2) 
El 2=2 .0D0*K1 13*Z ( 1 ) *Z (3) +K133* (Z (3) **2) 

QJAC (3 , 1 ) =E1 1+E1 2+PK4V 

Equation 3.3i 

E21 =2 .0D0*K1 22*Z ( 1 ) *Z (2) + (K2+M*ZG*UND) +3 .0D0*K222* (Z(2) **2) 
E22=2.0D0*K223*Z(2)*Z(3)+K233*(Z(3)**2)+K112*(Z(1 )**2) 

QJAC (3 ,2)=E21+E22+PK4R 

Equation 3.3k 

DGZ1 =GM*DC0S(Z (3) ) +0 . 5D0*BM* (DTAN (Z(3) ) ) **2*DC0S(Z (3) ) 
DGZ2=BM*(DTAN(Z(3) )*( (1 .0D0/DC0S(Z(3) ) )**2)*DSIN(Z(3) ) ) 
DGZ=DGZ1 +DGZ2 
DGZ=GM*DC0S(Z(3) ) 

Equation 3.3) 
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E31 =K3+K1 13* (Z( 1 ) **2) +2 .0D0*K133*Z ( 1 ) *Z(3) +K223* (Z(2) **2) 

E32=2 .0D0*K233*Z(2) *Z(3) -WG*DGZ 

* 

QJAC (3 , 3) =E31 +E32 

* 

RETURN 

END 

★ 

SUBROUTINE EIGTH1 (VDEL , V , R , RPHI ) 

★ 

INCLUDE ’ HYDCOMN.for’ 

INTEGER I , J , IERR 

REAL*8 AE(4,4) ,BE(4,4) ,BETA(4) ,VDEL(300) ,V(300) ,R(300) ,RPHI(300) 
REAL*8 P(300) , EVALR (4) , EVALI (4 ) , EVAL (4 ) ,ALFR(4) ,ALFI(4) ,Z(4,4) 
EXTERNAL DMACH , RGG 
DO 500 1=1 , IMAX 
P( I ) =0 . ODO 

* SET UP THE B MATRIX AND RESET EACH I 

UP=UND* (WP+TAU* ( ( V( I ) +XP*R ( I ) ) **2+CPV*V ( I ) +CPR*R( I ) ) ) 

V9=V ( I ) *0 .99D0 
R9=R ( I ) *0 . 99D0 

UPV9=UND* (WP+TAU* ( (V9+XP*R ( I ) ) **2+CPV*V9+CPR*R ( I ) ) ) 

UPR9=UND* (WP+TAU* ( (V( I ) +XP*R9) **2+CPV*V( I ) +CPR*R9) ) 

JP=(U*UP) / (NP*D) 

JPV9=(U*UPV9)/(NP*D) 

JPR9=(U*UPR9)/(NP*D) 

KT=KTA+(KTB*JP) 

KTV9=KTA+ (KTB*JPV9) 

KTR9=KTA+(KTB*JPR9) 

UR=UP*EPSILON*( (1 .0D0+8.0D0*K*KT)/(PI*(JP**2) ) )**0.5 
URV9=UPV9* EPSILON* ( ( 1 .0D0+8 . 0D0*K*KTV9) / ( PI* ( JPV9**2) ) ) **0 . 5 
URR9=UPR9* EPSILON* ( ( 1 .0D0+8 . 0D0*K*KTR9) / ( PI* ( JPR9**2) ) ) **0 . 5 
VR=GAMMA*V ( I )+CDR*R ( I )+CDRRR* (R ( I ) **3)+CDRRV*V( I ) * (R ( I ) **2) 
VRV9=GAMMA*V9+CDR*R ( I ) +CDRRR* (R ( I ) **3 ) +CDRRV*V9* ( R ( I ) **2 ) 
VRR9=GAMMA*V( I ) +CDR*R9+CDRRR* (R9**3) +CDRRV*V( I ) * (R9**2) 
VD=(UR**2+VR**2)**0.5 
VDV9=(URV9**2+VRV9**2) **0.5 
VDR9=(URR9**2+VRR9**2) **0.5 
ALPHAR=VDEL( I )+DATAN(VR/UR) 

ALPHARV=VDEL ( I ) + DAT AN (VRV9/URV9 ) 

ALPHARR=VDEL( I )+DATAN (VRR9/URR9) 

FLAM= (6.1 3D0* LAMBDA ) / ( LAMBDA+2 . 25D0 ) 

FNN=FLAM* ( AR / (L**2) ) * ( (VD) **2) *DSIN(ALPHAR) 

FNNV9=FLAM*(AR/(L**2) )*( (VDV9) **2) *DSIN(ALPHARV) 

FNNR9=FLAM* (AR/ (L**2) ) * ( (VDR9) **2) *DSIN(ALPHARR) 

★ 

Y4=- ( 1 .ODO+AH) *FNN*DCOS(VDEL( I ) ) 

Y4V9=- (1 .ODO+AH) *FNNV9*DC0S(VDEL( I ) ) 

Y4R9=- (1 .ODO+AH) *FNNR9*DC0S(VDEL ( I ) ) 

PY4V= (Y4 - Y4V9) / (V ( I ) - V9 ) 

PY4R=(Y4-Y4R9) / (R ( I ) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS(VDEL ( I ) ) 

N4V9=- (XR+AH*XH) *FNNV9*DC0S(VDEL ( I ) ) 

N4R9=- (XR+AH*XH) *FNNR9*DC0S(VDEL ( I ) ) 

PN4V= (N4-N4V9) / (V( I ) -V9) 

PN4R=(N4-N4R9) / (R ( I ) -R9) 

K4= ( 1 .ODO+AH) *ZR*FNN*DC0S( VDEL ( I ) ) 

K4V9= ( 1 .ODO+AH) *ZR*FNNV9*DC0S( VDEL ( I ) ) 

K4R9=( 1 .ODO+AH) *ZR*FNNR9*DC0S(VDEL ( I ) ) 

PK4V= (K4-K4V9) / (V( I ) - V9) 
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PK4R= (K4 -K4R9) / (R ( I ) -R9) 

Calculate the fundamental circular frequency 

WN=( (WG*GM) / (0.000021 DO) ) **0.5 
IF(FN.LE.0.1 )THEN 
KAPPA=0 . 1 DO 

ELSE I F ( FN . GE . 0 . 2 ) THEN 

KAPPA=0 . 2D0 

ELSE 

KAPPA=FN 

ENDIF 

KPHID0T=-0 . 000021 DO*KAPPA*WN 

Equation 4.1b 

BE ( 1 , 1 ) =0 . 00035245D0 

BE ( 1 , 2)=0 .000875D0 

BE(1 ,3)=0.0D0 

BE ( 1 ,4)=-0 .000213D0 

BE (2 , 1 ) =0 .014969D0 

BE ( 2 , 2 ) =0 . 00035245D0 

BE(2,3) = -0 .000221 DO 

BE(2,4)=0.0D0 

BE (3 , 1 ) = -0 .000221 DO 

BE(3,2)=0.0D0 

BE ( 3 , 3 ) =0 . 000021 DO 

BE ( 3 , 4 ) = - KPH IDOT 

BE (4 , 1 ) =0 .ODO 

BE(4,2)=0.0D0 

BE (4 ,3) =0 .ODO 

BE (4 ,4 )=1 .ODO 

Equation 4.1c SET UP THE A MATRIX 

C11=N1+3*N111*(V(I)**2)+N122*(R(I) **2) +2 .0D0*N1 1 2*V ( I ) *R ( I ) 

Cl 2=2 .0D0*N1 13*V( I ) *RPHI ( I )+N133* (RPHI ( I ) * *2 ) 

AE ( 1 , 1 ) =C1 1 +C1 2+PN4V 

C21 =2 .0D0*N122*V ( I ) *R ( I ) + (N2-M*XG*UND) +3 .0D0*N222* ( R ( I ) **2) 
C22=N1 12* (V ( I ) **2) +2 . 0D0*N223*R ( I ) * R PH 1(1) +N233* (RPHI ( I ) * * 2 ) 

AE(1 ,2)=C21+C22+PN4R 
AE ( 1 ,3)=0 .ODO 

C41=N3+N1 13*(V( I) **2) +2 .0D0*N133*V ( I ) * R PH 1(1) +N223* (R ( I ) **2) 
C42=2 .0D0*N233*R ( I ) *RPHI ( I ) 

AE ( 1 ,4)=C41+C42 

D1 1 =V1 +3 .0D0*Y1 11*(V(I)**2)+Y122*(R(I)**2)+2. 0D0*Y1 1 2*V ( I ) *R ( I ) 
D1 2=2 .ODO* Y1 13*V( I ) *RPHI ( I ) +Y133* (RPHI ( I ) * *2 ) 

AE (2 , 1 ) =D1 1 +D12+PY4V 

D21 =2 .0D0*Y1 22*V ( I ) *R ( I ) + (Y2-M*UND) +3 .0D0*Y222* (R ( I ) **2) 

D22=2 .0D0*Y223*R ( I ) *RPHI ( I ) +Y233* (RPHI ( I ) **2) +Y1 1 2* (V ( I ) **2 ) 

AE(2,2)=D21+D22+PY4R 
AE (2 ,3) =0 .ODO 



o o o 



* 



D41=Y3+Y1 13*(V(I)**2)+2 .0D0*Y133*V ( I ) *RPHI ( I ) +Y223* ( R ( I ) **2) 

D42=2 .0D0*Y233*R ( I ) *RPHI ( I ) 

AE(2,4)=D41+D42 

El 1 =K1 +3 .0D0*K1 11*(V(I)**2)+K122*(R(I)**2) +2 .0D0*K1 1 2*V ( I ) *R ( I ) 
E12=2 .0D0*K1 13*V(I)*RPHI(I) +K133* (RPHI ( I ) **2) 

AE (3 , 1 ) =E1 1 +E1 2+PK4V 

E21 =2 . 0D0*K1 22*V( I ) *R ( I ) + (K2+M*ZG*UND) +3 .0D0*K222* (R(I)**2) 

E22=2 . 0D0*K223*R ( I ) * R PH 1(1) +K233* (RPHI ( I ) **2) +K1 12*(V(I)**2) 

AE (3 , 2) =E21 +E22+PK4R 
AE(3,3)=0.0D0 

DGZ1 =GM*DCOS(RPHI ( I ) ) +0 .5D0*BM* (DTAN (RPHI ( I ) ) ) **2*DC0S(RPHI ( I ) ) 
DGZ2=BM*(DTAN(RPHI(I) )*( (1 .0D0/DC0S(RPHI ( I ) ) ) **2) *DSIN (RPHI ( I ) ) ) 
DGZ=DGZ1 +DGZ2 
DGZ=GM*DCOS(RPHI ( I ) ) 

E41 =K3+K1 13* ( V( I ) **2 ) +2 .0D0*K133*V ( I ) *RPHI ( I ) +K223* (R ( I ) **2) 

E42=2 . 0D0*K233*R ( I ) *RPHI ( I ) - WG*DGZ 

AE (3 , 4) =E41 +E42 
AE(4 , 1 ) =0 . 0D0 
AE(4,2)=0.0D0 
AE(4,3)=1 .ODO 
AE(4,4)=0.0D0 

* This is a subroutine that utilizes the EISPACK LIBRARY [Ref. 20) 



Purpose: Calls the recommended sequence of subroutines from the 

Eigensystem subroutine package (EISPACK) to find the 
eigenvalues of the real generalized eigenproblem : 

AX = ( LAMBDA) *BX 

Arguments : 

NM -Row dimension of the two-dimensional array parameters as 
declared in the calling program Dimension statement 
( Input) 

N -the order of the matrices A and B. (Input) 

A -Contains a real general matrix (Input) 

B -Contains a real general matrix (Input) 

MATZ -Path indicator. If= 0, eigenvalues are found. (Input) 
ALFR -Real parts of the numerators of the eigenvalues. (Output) 
ALFI -Imaginary parts of the numerators of the eigenvalues. 
(Output) 

BETA -Contains the denominators of the eigenvalues, which are 
given by the ratios: (ALFR+I*ALFI ) /BETA. (Output) 

Z -For eigenvector usage if MATZ not equal to zero. (Output) 
IERR -Error completion code. (Output) 



CALL RGG(4, 4, AE, BE, ALFR, ALFI, BETA, 0,Z, IERR) 
DO 600 J=1 ,4 
IF(BETA(J) .NE.0.0D0)THEN 
EVALR(J) =ALFR ( J ) /BETA( J ) 

EVALI ( J ) =ALFI ( J ) /BETA( J ) 
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ELSE 

EVALR ( J ) =DMACH ( 2 ) 

EVAL I ( J ) =DMACH ( 2 ) 

ENDIF 

600 CONTINUE 

WR ITE ( 25 , 25 ) EVALR ( 1 ) , EVAL 1(1), EVALR ( 2 ), EVAL 1(2) 

WRITE (30, 25) EVALR (3) , EVAL I (3) , EVALR (4 ) , EVAL I (4 ) 

* 

* Calculate the damping ratio and the natural frequency 

★ 

Z1 =- EVALR ( 1 ) / (DSQRT (EVAL I ( 1 ) **2+EVALR ( 1 ) **2 ) ) 

Z2=- EVALR (2) / (DSQRT ( EVAL I (2) **2+EVALR (2) **2) ) 

Z3=- EVALR (3) / (DSQRT (EVAL I (3) **2+EVALR (3) **2) ) 

Z4=- EVALR (4 ) / (DSQRT (EVAL I (4 ) **2+EVALR (4) **2) ) 

WN1 =DSQRT ( EVAL I ( 1 ) **2+EVALR ( 1 ) * * 2 ) 

WN2=DSQRT ( EVAL 1(2) * *2+EVALR ( 2 ) * * 2 ) 

WN3=DSQRT ( EVAL 1(3) * *2+EVALR ( 3 ) * *2 ) 

WN4=DSQRT ( EVALI (4) **2+EVALR ( 4 ) **2) 

WRITE (35 , 25) Z1 ,Z2,Z3,Z4 
WRITE(36,25)WN1 ,WN2,WN3,WN4 
25 FORMAT (4E15.5) 

500 CONTINUE 
RETURN 
END 

★ 

SUBROUTINE EIGTH2(VDEL , V , R ) 

* 

INCLUDE ’ HYDCOMN . f or ’ 

INTEGER I , J , IERR 

REAL*8 AE(2,2) ,BE(2,2) ,BETA(2) ,VDEL(300) ,V(300) ,R(300) 

REAL*8 EVALR (2) , EVAL I (2) , EVAL (2) , ALFR (2) ,ALFI (2) ,Z(2,2) 
EXTERNAL DMACH , RGG 
DO 500 I=1,imax 

* SET UP THE B MATRIX AND RESET EACH I 

UP=UND* (WP+TAU* ( (V ( I ) +XP*R ( I ) ) **2+CPV*V ( I ) +CPR*R ( I ) ) ) 

V9=V ( I ) *0 . 99D0 
R9=R ( I ) *0 . 99D0 

UPV9=UND* (WP+TAU* ( ( V9+XP*R ( I ) ) **2+CPV*V9+CPR*R ( I ) ) ) 

UPR9=UND* (WP+TAU* ( (V ( I ) +XP*R9) **2+CPV*V ( I ) +CPR*R9 ) ) 
JP=(U*UP)/(NP*D) 

JPV9=(U*UPV9) / (NP*D) 

JPR9= (U*UPR9) / (NP*D) 

KT=KTA+ (KTB*JP) 

KTV9=KTA+ ( KTB* JPV9 ) 

KTR9=KTA+ ( KTB* J PR9 ) 

UR=UP*EPSILON* ( (1 .0D0+8 .ODO*K*KT) / (PI*(JP**2) ) ) **0.5 
URV9=UPV9*EPSIL0N* ( ( 1 .0D0+8.0D0*K*KTV9) / ( PI* ( JPV9**2) ) ) **0 . 5 
URR9=UPR9*EPSIL0N* ( ( 1 .0D0+8 .0D0*K*KTR9) / ( PI* ( JPR9**2) ) ) **0 .5 
VR=GAMMA*V ( I ) +CDR*R ( I ) +CDRRR* ( R ( I ) * *3 ) +CDRRV*V ( I) * (R ( I ) **2) 
VRV9=GAMMA*V9+CDR*R ( I ) +CDRRR* (R ( I ) **3) +CDRRV*V9* (R ( I ) **2) 
VRR9=GAMMA*V ( I ) +CDR*R9+CDRRR* (R9**3) +CDRRV*V ( I ) * (R9**2) 

VD= (UR**2+VR**2) **0 .5 
VDV9= (URV9**2+VRV9**2) **0.5 
VDR9= (URR9**2+VRR9**2) **0.5 
ALPHAR=VDEL ( I ) +DATAN ( VR / UR ) 

ALPHARV=VDEL ( I ) +DATAN (VRV9/URV9 ) 

ALPHARR=VDEL ( I ) +DATAN (VRR9/URR9) 

FLAM= (6.1 3D0* LAMBDA ) / ( LAMBDA+2 . 25D0 ) 

FNN=FLAM*(AR/(L**2) )*( (VD) **2) *DSIN (ALPHAR) 

FNNV9=FLAM* (AR/(L**2))*( ( VDV9) **2) *DSIN (ALPHARV) 

FNNR9=FLAM* ( AR/ (L**2) ) * ( ( VDR9) **2) *DSIN (ALPHARR ) 
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Y4=- ( 1 .ODO+AH) *FNN*DCOS(VDEL( I) ) 

Y4V9=- ( 1 . ODO+AH) *FNNV9*DC0S(VDEL ( I ) ) 

Y4R9=- (1 .ODO+AH) *FNNR9*DC0S(VDEL ( I ) ) 

PY4V= ( Y4 - Y4V9 ) / (V ( I ) -V9) 

PY4R= (Y4- Y4R9) / ( R ( I ) -R9) 

N4=- (XR+AH*XH) *FNN*DCOS(VDEL ( I ) ) 

N4V9=- (XR+AH*XH) *FNNV9*DC0S(VDEL( I) ) 

N4R9=- (XR+AH*XH)*FNNR9*DC0S(VDEL( I) ) 

PN4V=(N4-N4V9) / (V(I) -V9) 

PN4R=(N4-N4R9) / (R(I) -R9) 

WN=( (WG*GM)/ (0.000021 DO) ) * *0.5 

Equation 4.4b 

BE ( 1 , 1 ) =0 .00035245D0 
BE( 1 ,2)=0 .000875D0 
BE ( 2 , 1 ) =0 . 0 1 4969D0 
BE (2 , 2) =0 .00035245D0 

Equation 4.4c-g 

C11=N1+3*N111*(V(I)**2)+N122*(R(I)**2)+2.0D0*N1 12*V(I)*R(I) 

AE ( 1 , 1 ) =C1 1 +PN4V 

C21 =2 .0D0*N1 22*V ( I ) *R ( I ) + (N2-M*XG*UND) +3 .0D0*N222* (R ( I ) **2) 
C22=N1 1 2* (V( I ) **2) +PN4R 

AE ( 1 ,2)=C21+C22 

Dll =Y1 +3 .0D0*Y1 11*(V(I) A *2)+Y122*(R(I)**2)+2 .0D0*Y1 1 2*V ( I ) *R ( I ) 
AE (2 , 1 ) =D1 1 +PY4V 

D21 =2 .0D0*Y1 22*V( I ) *R ( I ) + (Y2-M*UND) +3 .0D0*Y222* (R ( I ) **2) 

D22=Y1 12* (V( I ) **2) +PY4R 

AE(2,2)=D21+D22 

CALL RGG(2,2 > AE > BE,ALFR > ALFI > BETA > 0 > Z I I ERR) 

DO 600 J=1 ,2 

IF(BETA(J) .NE.O .ODO)THEN 
EVALR(J)=ALFR(J) /BETA(J) 

EVALI(J)=ALFI(J) /BETA(J) 

ELSE 

EVALR ( J ) =DMACH ( 2 ) 

EVAL I ( J ) =DMACH ( 2 ) 

ENDIF 

600 CONTINUE 

WRITE (26 , 25) EVALR ( 1 ) , EVALI ( 1 ) , EVALR (2) , EVAL I (2) 

25 FORMAT (4E15.5) 

500 CONTINUE 
RETURN 
END 

SUBROUTINE EIGTH3(VDEL, RPHI ) 

INCLUDE ’ HYDCOMN.for’ 

INTEGER I , J , IERR 
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REAL*8 AE(2,2) ,BE(2,2) ,BETA(2) ,VDEL(300) ,RPHI(300) 

REAL*8 P (300) ,EVALR(2) ,EVALI(2) ,EVAL(2) ,ALFR(2) ,ALFI(2) , 1 ( 2 , 2 ) 
EXTERNAL DMACH, RGG 
DO 500 I=1,imax 
P( I)=O.ODO 

* SET UP THE B MATRIX AND RESET EACH I 

WN=( (WG*GM) / (0.000021 DO) ) **0.5 
IF(FN.LE.O. 1 )THEN 
KAPPA=0 . 1 DO 

ELSE IF(FN.GE.0.2)THEN 

KAPPA=0 . 2D0 

ELSE 

KAPPA=FN 

ENDIF 

KPHID0T=-0.000021D0*KAPPA*WN 

* 

* Equation 4.7b 

* 

BE(1 ,1 )=0.000021D0 
BE( 1 ,2)=-KPHID0T 
BE (2 , 1 ) =0 .OdO 
BE(2,2)=1 .ODO 

★ 

* Equation 4.7c-e 

★ 

AE(1 ,1 )=O.OdO 

c DGZ1=GM*DC0S(RPHI ( I ) ) +0 . 5D0*BM* (DTAN (RPHI ( I ) ) ) **2*DC0S(RPHI ( I ) ) 

c DGZ2=BM* (DTAN (RPHI ( I ) ) * ( ( 1 .ODO/DCOS(RPHI ( I) ) ) **2) *DSIN (RPHI ( I ) ) ) 

c DGZ=DGZ1 +DGZ2 

DGZ=GM*DCOS(RPHI ( I) ) 

AE ( 1 ,2)=K3-WG*DGZ 
AE(2,1 )=1 .ODO 
AE(2,2)=0.0D0 

CALL RGG(2,2,AE,BE,ALFR,ALFI ,BETA,0,Z, I ERR) 

DO 600 J=1 ,2 
IF(BETA(J) .NE.0.0D0)THEN 
EVALR ( J ) =AL FR ( J ) / BETA ( J ) 

EVALI (J)=ALFI (J) /BETA(J) 

ELSE 

EVALR ( J ) =DMACH ( 2 ) 

EVAL I ( J ) =DMACH ( 2 ) 

ENDIF 

600 CONTINUE 

WRITE (27,25) EVALR ( 1 ) , EVALI ( 1 ) , EVALR (2) ,EVALI(2) 

Z1 = - EVALR ( 1 ) / ( DSQRT ( EVAL I ( 1 ) * * 2+ EVALR ( 1 ) * *2 ) ) 

Z2= - EVALR ( 2 ) / ( DSQRT ( EVAL I ( 2 ) * * 2+ EVALR ( 2 ) * *2 ) ) 

WN1 =DSQRT ( EVALI ( 1 ) * *2+EVALR ( 1 ) * *2 ) 

WN2=DSQRT (EVALI (2) **2+EVALR(2) **2) 

WRITE(37,25)Z1 ,Z2,WN1 ,WN2 
25 FORMAT (4E1 5 . 5) 

500 CONTINUE 
RETURN 
END 
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APPENDIX B 



Fn = 0.3 and a — 0.6 
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Fn = 03 m and a = 0.6 




B.2: Effect of GM on Coupled Model Roll Eigenvalues; Imaginary Part 
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GM = 0.3 m and a = 0.6 




B.3: Effect of GM on Coupled Model Upper Steering Eigenvalues 



85 



GM = 0.3 m and a = 0.6 
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B.4: Effect of GM on Coupled Model Lower Steering Eigenvalues 
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GM = 0.3 m and a - 0.6 
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GM = 0.3 m and Fn = 0.6 
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GM = 0.3 in and Fn = 0.6 
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B.7: Effect of Fn on Coupled Model Upper Steering Eigenvalues 
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GM = 0.3 m and a = 0.6 
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B.8: Effect of Fn on Coupled Model Lower Steering Eigenvalues 
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RUDDER ANGLE (deg) 



GM = 0.3 m and a = 0.6 




B.9: Effect of Fn on Coupled Model Roll Eigenvalues 
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GM = 0.3 m and a = 0.6 
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3.10: Effect of Fn on Coupled Model Root Locus 
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GM = 0.3 m and Fn = 0.3 




B.ll: Effect of Speed Loss on Coupled Model Roll Eigenvalues 
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GM = 0.3 m and Fn = 0.3 




3.12: Effect of Speed Loss on Coupled Model Root Locus 
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GM = 0.3 m and Fn = 0.3 
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3.13: Effect of Speed Loss on Coupled Model Upper Steering 

Eigenvalues 
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GM = 0.3 m and Fn = 0.3 
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B.14: Effect of Speed Loss on Coupled Model Lower Steering 

Eigenvalues 
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Fn = 0.3 
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B.15: Effect of Speed Loss on Coupled Model Sway Velocity with GM = 1.0 m 



97 



RUDDER ANGLE (deg) 



Fn = 0.3 




A1I3013A AVMS * 



B.16: Effect of Speed Loss on Coupled Model Sway Velocity with GM = 3.0 m 
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Fn = 0.3 
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B.17: Effect of Speed Loss on Coupled Model Roll Angle with GM = 1.0 m 
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B.18: Effect of Speed Loss od Coupled Model Roll Angle with GM = 3.0 m 
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Fn = 0.3 
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B.19: Effect of Speed Loss on Coupled Model Root Locus with GM = 1.0 m 
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B.20: Effect of Speed Loss on Coupled Model Root Locus with GM = 3.0 m 
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